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Preface 


This book - an in-depth examination of chemical thermodynamics - is 
written for an audience of engineering undergraduates and Masters students 
in the disciplines of chemistry, physical chemistry, process engineering, 
materials, etc., and doctoral candidates in those disciplines. It will also be 
useful for researchers at fundamental- or applied-research labs, dealing with 
issues in thermodynamics during the course of their work. 

These audiences will, during their undergraduate degree, have received a 
grounding in general thermodynamics and chemical thermodynamics, which 
all science students are normally taught, and will therefore be familiar with 
the fundamentals, such as the principles and the basic functions of 
thermodynamics, and the handling of phase- and chemical equilibrium 
states, essentially in an ideal medium, usually for fluid phases, in the absence 
of electrical fields and independently of any surface effects. 

This set of books, which is positioned somewhere between an 
introduction to the subject and a research paper, offers a detailed 
examination of chemical thermodynamics that is necessary in the various 
disciplines relating to chemical- or material sciences. It lays the groundwork 
necessary for students to go and read specialized publications in their 
different areas. It constitutes a series of reference books that touch on all of 
the concepts and methods. It discusses both scales of modeling: microscopic 
(by statistical thermodynamics) and macroscopic, and illustrates the link 
between them at every step. These models are then used in the study of solid, 
liquid and gaseous phases, either of pure substances or comprising several 
components. 
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The various volumes of the set will deal with the following topics: 

- phase modeling tools: application to gases; 

- modeling of liquid phases; 

- modeling of solid phases; 

- chemical equilibrium states; 

- phase transformations; 

- electrolytes and electrochemical thermodynamics; 

- thermodynamics of surfaces, capillary systems and phases of small 
dimensions. 

Appendices in each volume give an introduction to the general methods 
used in the text, and offer additional mathematical tools and some data. 

This series owes a great deal to the feedback, comments and questions 
from all my students are the Ecole nationale superieure des mines 
(engineering school) in Saint Etienne who have “endured” my lecturing in 
thermodynamics for many years. I am very grateful to them, and also tha nk 
them for their stimulating attitude. This work is also the fruit of numerous 
discussions with colleagues who teach thermodynamics in the largest 
establishments - particularly in the context of the group “Thermodic”, 
founded by Marc Onillion. My thanks go to all of them for their 
contributions and conviviality. 

This volume in the series is devoted to the study of liquid phases. 

Chapter 1 describes the modeling of pure liquids, either using the radial 
distribution function or partition functions. The different models presented 
herein range from the very simplest to the most complex. The results yielded 
by these models are then compared, both to one another and to the results 
found experimentally. 

The second chapter describes the tools used for macroscopic modeling of 
solutions. The use of limited expansions of the activity coefficient logarithm 
is presented, before we define simple solution models such as the ideal dilute 
solution, regular solutions and athermal solutions, on the basis of 
macroscopic properties. 
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Next, in Chapter 3, we present a number of solution models with 
microscopic definition, including random distribution models and models 
integrating the concepts of local composition and combinatorial excess 
entropy. 

The fourth chapter deals with the modeling of ionic solutions combining 
the term due to the electrical effects, found using the Debye and Hiickel 
model, with the terms of local composition and combinatorial excess entropy 
found in the previous chapter. 

Chapter 5 presents the various experimental methods for determining the 
activity or the activity coefficient of a given component in a solution. 

Finally, three appendices are provided, which recap a few notions about 
statistical methods of numerical simulation (Appendix 1), and offer some 
reminders about the properties of solutions (Appendix 2) and statistical 
thermodynamics (Appendix 3) - subjects which were discussed in detail in 
the first volume of this series. 


Michel SOUSTELLE 
Saint- Vallier, 
April 2015 




Notations and Symbols 


{gas} pure, {{gas}} in a mixture, (liquid) pure, ((liquid)) in solution, (solid) pure, 
((solid)) in solution 

A: area of a surface or an interface. 

(12) 

d|| : Hamaker constant between two media 1 and 2. 

A: affinity. 

A : electrochemical affinity. 

A m '. molar area. 

A,„: molecular area. 

a: cohesion pressure of a gas or radius of the unit cell of a liquid. 

A, B, . . . : components of a mixture. 

a"" x and b m,x \ mixing terms of the constants in a state equation. 

B : 1 th coefficient of the virial in the pressure expansion. 

Bf. 1 th coefficient of the virial. 

b\ covolume of a gas or cosurface of an adsorbed gas. 

C: concentration or concentration in a potential-pH plot. 

nxs . 

• 


molar heat capacity of excess at constant pressure. 
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Q: 

C+ : 

Cv(el)- 

C v ( r ). 


c, 


v(t)- 


c 


v(v)- 


Cv, Cp. 


D : 

D(T/0 d ): 

d: 

A e S: 

d,: 

djS: 

dco: 

E: 

E: 

E(T/0 e ): 

E 0 : 

£°: 

-^abs- 

Eb: 


molar concentration (or molarity) of component i. 

mean concentration of ions in an ionic solution. 

contribution of free electrons in a metal to the molar heat capacity. 

contribution of rotational motions to the heat capacity at constant 
volume. 

contribution of translational motions to the heat capacity at 
constant volume. 

contribution of vibrational motions to the heat capacity at constant 
volume. 

heat capacity at constant volume and constant pressure, 
respectively. 

capacity of a capacitor or number of independent components. 

dielectric constant of the medium or diameter of protection or 
contact of a molecule. 

Debye’s function. 

distance between two liquid molecules, 
entropy exchange with the outside environment, 
degree of oxidation i of an element A. 
internal entropy production, 
elementary volume, 
energy of the system. 

Young’s modulus. 

Einstein’s function. 

internal energy associated with a reaction at a temperature of OK. 

standard electrical potential or standard electromotive force (emf) 
of an electrochemical cell. 

reversible emf of an electrochemical cell. 

balance equation. 
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E c : 


E,: 

E f . 

Ekin' 

E„: 


e 0 : 


^abs- 


F: 

riinix . 

"m 


— mix 
Fi 


Ft: 

F el : 

F : 

F m : 

/het • 

ft- 

fi- 

for f?: 


mean total energy of an element in the canonical ensemble. 

total energy of the canonical ensemble. 

potential energy due to interactions. 

energy of an element j of the canonical ensemble. 

molar kinetic energy of electrons in a metal. 

set of variables with p intensive variables chosen to define a 
system. 

relative emf of an electrode, 
standard emf of an electrode. 

equi-activity- or equiconcentration emf of an electrode, 
absolute emf of an electrode. 

Helmholtz energy. 

molar excess Helmholtz energy. 

partial molar excess Helmholtz energy of the component i. 

partial molar mixing Helmholtz energy of the component i. 
free energy, partial molar Helmholtz energy of the component i. 
contribution of free electrons to the molar Helmholtz energy. 

electrochemical Helmholtz energy, 
molar Helmholtz energy, 
faraday. 

heterogeneous wetting function, 
fugacity of the component i in a gaseous mixture, 
molar Helmholtz energy of pure component i. 
fugacity of a pure gas i. 
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G ™ : excess Gibbs energy. 

G a : electrocapillary Gibbs energy. 

G : electrochemical Gibbs energy. 

— xs 

Gi : partial excess molar Gibbs energy of component i. 

G, Gi , [G]: free enthalpy, partial molar free enthalpy of i, generalized free 

enthalpy. 

G m : molar Gibbs energy. 


G mix . 
m 

g- 
g!’ ■ 

ga- 

gi- 

g(e> 

g( r): 

g(vj- 

g- 

Hj : 


molar Gibbs energy of mixing. 

osmotic coefficient or acceleration due to gravity or degeneration 
coefficient or multiplicity or statistical mass. 

molar Gibbs energy of pure component i. 

statistical weight of fundamental electron level of nucleus a. 

coefficient of multiplicity of state i. 

statistical weight of electron levels. 

radial distribution function. 

distribution of velocity components along Ox axis. 

molar Gibbs energy of gas i at pressure of 1 atmosphere in a 
mixture. 

standard molar enthalpy of formation at temperature T. 


H , Hi : enthalpy, partial molar enthalpy of i. 

H: Hamiltonian. 

Hj\ integral of resonance between two neighboring identical atoms. 

H u : Coulombian integral between two neighboring identical atoms. 

5C. magnetic field. 

fj '■ electrochemical enthalpy. 
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H™ : molar excess enthalpy. 

H™ lx : molar mixing enthalpy. 

— XS 

Hi : partial excess molar enthalpy of component i. 

— mix 

Hi : partial molar mixing enthalpy of component i. 

h, : spreading coefficient. 

h: stoichiometric coefficient of protons in an electrochemical 

reaction. 

h: Planck’s constant. 

hf : molar enthalpy of pure component i. 

h sp \ Harkins spreading coefficient of a liquid on another. 

/: ionic strength of a solution of ions. 

I m : ionic strength in relation to molality values. 

/, / /, /? I 3 \ moments of inertia. 

If. integral of configuration of the canonical distribution function of 

translation. 

Van’t Hoff factor. 

Jj : partial molar value of/ relative to component i. 

J'j' ux : mixing value of/ relative to component i. 

j" nx : partial molar mixing value of /relative to component i. 

* 

Jj : value of /relative to component i in a perfect solution. 

its 

Jj : partial molar value of/ relative to component i in a perfect 

solution. 

fj : value of J for the pure component i in the same state of 

segregation. 
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K hj ^) : 

: 


K:, 


-^ads- 

Kax- 

A a P) . 

A i 


K d : 

K k : 

K { r c) : 

KP: 

K IP) ■ 

K r : 

K s \ 

k: 


k B : 


L t : 


lc. 

M: 

M: 

T: 

nr. 


rotational quantum number. 

thermodynamic coefficient associated with the set of variables E/i. 
Xj is its definition variable and Y t its definition function. 

constant of change of equilibrium for phase transition Tr for 
component i. 

weighting factor of local composition. 

equilibrium adsorption constant, 
solubility product of solid AX. 

coefficient of sharing of compound i between the two phases a 
and (3. 

dissociation constant, 
adsorption equilibrium function. 

equilibrium constant relative to concentrations, 
equilibrium constant relative to fugacity values. 

equilibrium constant relative to partial pressure values, 
equilibrium constant, 
solubility product, 
wavenumber. 

Boltzmann’s constant. 

latent heat accompanying the transformation t. 
capillary length, 
molar mass. 

magnetic moment or Madelung constant, 
mass of solute s in grams per kg of solvent, 
total mass. 
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m{. 

N: 


N a : 


N c : 
N c : 


n 


(a). 


<n>\ 


N l : 


pinix . 

P: 

psubl . 
pvap pO . 


mass of component i. 

number of components of a solution or a mixture of gases or 
involved in a reaction or number of molecules of a collection. 

Avogadro’s number. 

number of molecules of component A. 

number of elements in the canonical collection. 

total number of cells of a liquid. 

number of objects i in the system with energy 6; or number of 
moles of component i. 

translational quantum nimber or total number of moles in a 
solution or a mixture. 

total number of moles in a phase a. 

mean number of neighboring vacancies of a molecule in a liquid, 
total number of vacancies in a liquid. 

critical pressure of the mixture, 
pressure of a gas. 

sublimating vapor pressure of component i. 
saturating vapor pressure of component i. 


P ' l '. nlx : relative pressure of the mixture. 

P c : critical pressure. 

Pf. partial pressure of component i. 

Pf. proportion of number of elements in a state j. 

p: number of external physical variables or spreading parameter. 

Fermi pulse, 
heat involved. 

reaction quotient in terms of activity. 


Pf- 

Q- 

Qa- 
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heat of transformation at constant pressure; quotient of reaction in 
terms of partial pressures. 

reaction quotient of transformation r. 

transformation heat at constant volume. 

equilibrium heat of adsorption. 

differential heat of adsorption. 

volumetric fraction parameter. 

isosteric heat of adsorption. 

reaction rate 

perfect gas constant. 

mean curvature radius of a surface or rate of reflux of distillation. 

radius of the ionic atmosphere. 

minimum distance of energy between two molecules. 

radius of a cylindrical tube. 

volumetric fraction parameter. 

Kelvin radius. 


molar mixing entropy. 


partial excess molar entropy of component i. 


partial mixing molar entropy of component i. 


oversaturation of a solution. 


entropy or partial molar entropy of i. 
electrochemical entropy. 


excess molar entropy. 


parameter of order of an alloy. 
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sf* : molar entropy of pure component i. 

T: temperature. 

T' c mx : critical temperature of the mixture. 



U, Ui : 
U el : 

U : 

U m : 

U R : 
u+, u 


second-order transition temperature, 
relative temperature of the mixture, 
boiling point of azeotropic solution, 
critical temperature. 

Fermi temperature, 
boiling point of pure i. 
melting point of pure i. 
sublimation temperature, 
vaporization temperature. 

excess molar internal energy, 
mixing molar internal energy, 
excess partial molar internal energy of component i. 

partial mixing molar internal energy of component i. 
internal energy, partial molar internal energy of i. 
contribution of free electrons to the molar internal energy, 
internal electrochemical energy. 

molar internal energy. 

crosslink internal energy. 

ionic mobilities of the cation and anion. 

molar internal energy of pure component i. 
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V, Vi : 

volume, partial molar volume of 

V c : 

critical volume. 

V G : 

Gibbs variance. 

V m : 

molar volume. 

V c \ 

volume of the unitary element of a liquid. 

v D - 

Duhem variance. 

V f- 

free volume per molecule. 

Voo: 

volume of influence around a molecule. 

vf: 

molar volume of pure component 

v: 

quantum vibration number. 

Vm- 

molecular volume. 

Vm- 

molar volume of solid at melting point. 

^mono • 

volume of monolayer of adsorbed gas. 

V x i- 

component along Ox axis of the velocity of a particle i. 

Wu- 

energy per square meter of interaction between the surfaces of 
phases 1 and 2. 

w t : 

mass fraction of the component i. 

w if . 

energy of exchange between atoms i and j. 

r («) . 
x k ■ 

molar fraction of component k in phase a. 

x,y,z: 

coordinates of a point in space. 

x,\ 

molar fraction of the component i in a solution. 

< y > - 

mean value of y. 

Yj and X{. 

conjugal intensive and extensive values. 

y<,/- 

Mayer function. 

yr- 

molar fraction of component i in a gaseous phase. 

Z: 

compressibility coefficient. 
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Z i- 
^ mix . 


Z c : 


z a a > : 

Z C(I)' 

Z C(I>- 

z: 


Z e - 

z;. 

z int- 


Z„\ 

Z pf- 

z,:. 


z{- 


z m- 


z v : 

or. 


or d . 

P 


rfe): 


compressibility coefficient of gas i. 

compressibility coefficient of the mixture of gases. 

molecular partition function of interaction between molecules. 

canonical partition function. 

canonical partition function of component A. 

canonical partition function of interaction. 

canonical partition function of translation. 

molecular partition function, altitude of a point or coordination 

index, number of nearest neighbors. 

electron molecular partition function or electrovalence of ion i. 

number of molecules that are near neighbors of a molecule i. 

contribution of internal motions to the molecular partition 
function. 

molecular partition function of nuclei, 
molecular partition function of a perfect gas. 
rotational molecular partition function, 
translational molecular partition function, 
translational molecular partition function of a perfect gas. 
vibrational molecular partition function. 

coefficient of dissociation of a weak electrolyte or linear dilation 
coefficient at pressure P or relative volatility or Lagrange 
multiplier relating to the number of objects of a collection or 
polarizability of a molecule. 

apparent dissociation coefficient of a weak electrolyte. 

Lagrange multiplier relating to the energy of the objects in a 
collection or volumetric dilation coefficient at pressure P. 

characteristic function with the set as canonical variables. 
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r k : coefficient of activity of a group. 

F. characteristic function. 

/]: surface excess or surface concentration of component i. 

r LJ : excess surface or surface concentration of component i in relation 

toy. 

y : coefficient of activity of the component i irrespective of the 

reference state or Griineisen parameter or structure coefficient 
whose value is x/2 for cubic crystal lattices with centered faces. 

Yo : activity coefficient of a solvent. 

Yf. activity coefficient of the species i or Griineisen factor of phonon 

i. 


rl UI) : 

r±- 

%■ 

Ao • 

4 - a t '■ 

A, A: 

Su- 

8 : 

£ A(A) ■ 


activity coefficient of component i, pure-substance reference. 

activity coefficient of component infinitely dilute solution 
reference. 

activity coefficient of component i, molar solution reference. 

mean activity coefficient of ions in an ionic solution, 
activity coefficient of a solute, 
spreading of a liquid. 

standard value at temperature T of A associated with the 
transformation r. 

value de A associated with the transformation r. 

Kronecker delta. 

coefficient of pressure increase at volume V. 
network energy of an atom of A in network A. 


Wagner interaction coefficient. 
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£ : 
£ 0 : 


£ c \ 


£/■ 


£/■'■ 


£ i(e>- 


£i(r)- 


£i(ty 


£i(v)- 

ZiS- 


^rnm* 

£o- 


£ p \ 

£ r ep' 

ij: 

%■ 

6b: 

6b: 

< 9 r : 

a 


electrical permittivity of the medium, 
electrical permittivity of a vacuum, 
energy of attraction between molecules, 
kinetic energy of a molecule, 
energy of the C-H bond. 

energy from the dispersion effect between molecules. 

Fermi energy. 

electronic energy of a molecule i. 
interactional energy of a molecule 
nuclear energy of a molecule i. 
rotational energy of a molecule i. 
translational energy of a molecule i. 
vibrational energy of a molecule i. 

energy of interaction between two molecules i and j or pair energy 
between atoms i and j. 

switch. 

energy due to the effect of orientation between molecules, 
potential energy of a molecule, 
repulsion energy between molecules, 
viscosity. 

Warren and Cowley’s order parameter. 

Debye’s vibration temperature. 

Einstein’s vibration temperature, 
characteristic rotation temperature, 
overlap fraction. 
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ec. 

k 
Ao+, 
A a ■ 
A, : 
A: 


A 0 : 

Mi, \Mk Mi ■ 




M- 

v. 


V H P ) '■ 


Vd- 


v e : 


$ 

n d : 

P- 

Pi 1 ')'- 

a 

o c : 

o*: 

T+, r_: 


surface fraction of a component, 
linear dilation coefficient. 

equivalent ionic conductivities of the cation and anion, 
absolute activity of component A. 
lateral chemical potential of component i. 

equivalent conductivity of an electrolyte or thermal wavelength of 
a molecule. 

maximum equivalent conductivity of an electrolyte. 

chemical potential of the component i, dipolar electrical moment 
of molecule i, generalized chemical potential. 

chemical potential of component i in liquid/gaseous state, 
respectively. 

electrochemical potential, 
vibration frequency. 

algebraic stoichiometric number of component A k in reaction p. 
Debye’s maximum frequency. 

stoichiometric coefficient of electrons in an electrochemical 
reaction. 

reaction extent. 

disjunction pressure. 

density of molecules in a spherical crown of radius r or volumetric 
density of electrical charges or density. 

density of molecules in an enclosure. 

surface energy or symmetry number. 

surface density of electrical charges. 

surface tension. 

cationic and anionic transport numbers. 
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0 : practical osmotic coefficient; expansion pressure. 

0j : coefficient of fugacity of component i in a gaseous mixture. 

0. coefficient of conductivity of a strong electrolyte or number of 

Phases. 

(pi : coefficient of fugacity of gas i in a mixture or volume fraction of a 

component. 

0° or 0 P : coefficient of fugacity of a pure gas. 

Xi : calorimetric coefficient relative to the variable x t . 

X'- electrical conductivity. 

Xt : coefficient of compressibility at temperature T. 

'Fj : electrostatic potential of ionic atmosphere. 

F(r) : electrostatic potential. 


F km : energy term between two groups, 

waveflinction. 

£2 be : number of complexions in Bose-Einstein statistics. 

Q . c : number of complexions in Fermi-Dirac statistics. 

il: number of complexions. 

at,: set of position coordinates of molecule 

OJx'. rotational velocity component in direction Ox. 
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Pure Liquids 


This chapter will be given over to atomic and molecular liquids. A pure 
molecular liquid is a liquid comprising only one type of non-dissociated 
molecules. The study of liquids is more difficult than that of gases and solids 
because they are in an intermediary state, structurally speaking. Indeed, as is 
the case with solids, we can imagine that in liquids (and this is confirmed by 
X-ray diffraction), the interactions between molecules are sufficiently 
powerful to impose a sort of order within a short distance of the molecules. 
However, the forces involved in these interactions are sufficiently weak for 
the molecules to have relative mobility and therefore for there to be disorder 
(no form of order) when they are far apart, as is the case with gases. 


1.1. Macroscopic modeling of liquids 

In the areas where liquids are typically used, far from the critical 
conditions, it is often possible to consider liquids to be incompressible - 
meaning that (OF / 3P) = 0 - but dilatable. The order of magnitude of a 

dilation coefficient is 10' 3 degrees’ 1 , whereas that of the compressibility 
coefficient is lO^atm’ 1 . 

As we approach the critical conditions, this approximation is no longer 
possible, and the properties of the liquid tend more to be governed by an 
equation of state. Whilst the “cubic” equations of state for gases do include 
critical conditions, it is accepted that the properties of liquids often 
necessitate equations of state that take account of the intervention of forces 
when more than two bodies are concerned. Additionally, the third- and 
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2 Modeling of Liquid Phases 


fourth coefficients of the virial, which can no longer be ignored in the case 
of liquids, become necessary when these types of forces are at work. 

Certain equations of state specific to liquids have been put forward in the 
literature, including Rocard’s, which is written thus: 


P = 


RTV 2 
(V-b/ 3) 3 


a 

v 2 


[ 1 . 1 ] 


In addition, this equation, expressed as the expansion of the virial, 
assumes the form: 


PV -RT 


1 + 


f h a ^ 1 

b 




RT 

ab 


b 

v 


ab b 


2 A 


RT 3 


V 2 


2 h 


27 3 RT 


1 ab 3 

— + - 


V 


1 


3 RT V 4 


[ 1 . 2 ] 


This equation does indeed include the third and fourth coefficients of the 
virial. 

The heat capacities at constant volume and constant pressure are 
practically identical, around 0.5cal/g, or2.1kJ/kg. 


1.2. Distribution of molecules in a liquid 

On a structural level, liquids are classified into two categories: associated 
liquids and non-associated liquids. 

A liquid is said to be non-associated if the intra-molecular degrees of 
freedom (rotational, vibrational, electronic and nuclear) are not majorly 
disturbed by the proximity of neighboring molecules. These liquids can be 
treated, as is the case with gases, with independence between the internal 
motions and the translation of the molecules. 

A liquid is said to be associated if, unlike in the previous case, the 
molecule’s internal degrees of freedom are disturbed by the proximity of 
other molecules. This disturbance may be so great that, in practical terms, we 
need to consider associations between molecules, coming together to form 
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dimers, trimers, etc. The new bonds that need to be taken into account are 
usually hydrogen bonds, whose energy is 4-5 times less than that involved in 
typical chemical bonds, but which are 4-5 times stronger than intermolecular 
bonding by van der Waals forces. When the temperature rises, these bonds 
are broken and, particularly when the thermal agitation energy (k B 7) is much 
greater than the energy in the hydrogen bond, the molecules separate and 
regain individuality when they are near to the gaseous state. 

These associations lend associated liquids very special properties, such as 
anomalies of the dilation coefficient, high viscosity, low surface tension and 
a high boiling point. Liquid water belongs to this category. The best way of 
dealing with these liquids in thermodynamics is to consider them no longer 
as pure liquids, but rather to treat them as associated solutions, with dimeric, 
trimeric (etc.) molecules - see section 2.5. 


1 . 2 . 1 . Molecular structure of a non-associated liquid 

Hereinafter, we shall focus only on non-associated liquids, and we shall 
suppose the molecules are spherical. A non-associated liquid is characterized 
by a local order, or short-distance order. The best illustration of this is of 
liquid metals. In Figure 1.1, which gives a 2-dimensional image of the 
arrangement of spherical molecules in a liquid, we can see that the molecules 
are relatively close together, and that around each molecule, there is an area 
of order which is illustrated by the circles superimposed on the figure. The 
short-distance arrangement, within the circles, is almost identical to the 
molecular arrangement in a solid crystal but, unlike with a crystal, there is no 
long-distance order. The two circles on Figure 1.1 exhibit no periodicity. 



Figure 1.1. Two-dimensional diagram of the distribution of molecules in a liquid 
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The second difference between a crystal and a liquid is that in the latter, 
the molecules are in perpetual motion, so Figure 1.1 is representative of the 
situation only at a given time. Unlike with a solid crystal, the distribution of 
those molecules would be different an instant later, although we would find 
similar zones of ordered arrangement. 

Flence, in order to accurately describe a liquid, we cannot content 
ourselves with merely describing the position of a few appropriately-chosen 
neighboring molecules, as we can with the lattice of the crystal. We would 
have to define the positions of each of the molecules at every moment in 
time. In view of the impossibility of the task in a medium with normal 
dimensions (around a mole, which contains 10 22 molecules), we use 
statistical methods using so-called correlation functions. The paired 
correlation function which we intend to examine constitutes the first level of 
this description. 


1 . 2 . 2 . The radial distribution function 

Throughout this chapter, we shall suppose that the interactions between N 
particles of a liquid medium are additive and paired, meaning that the 
internal energy due to these interactions is merely the sum of the interactions 
between molecules, two by two. Thus, the internal energy is the sum of the 
energies between the molecules taken two by two £,,(>',)• This energy 

depends only on the distance between the two molecules. Hence, we have: 

U(l,2,..JV) = fX;fc/) [1-3] 

i<j 


Consider a molecule chosen at random in the structure (Figure 1.2). Let 
d /V( /•) signify the number of molecules whose centers are situated in the 
crown between the two spheres centered on the chosen molecule, with radii r 
and r+dr and volume 471 r 2 dr. The density of molecules in the crown p(r), 
i.e. the number of molecules situated in the crown per unit volume of that 
spherical crown, at a distance r from the central molecule, is such that: 


f 

r 


1 d N(r) 
4 nr 2 dr 


p(r) = 


V, dV(r) J 


[1.4] 
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Figure 1.2. Arrangement of molecules of liquid around the center of a cage 


The volumetric density p is defined as the ratio of the total number of 
molecules in the liquid in question to the volume of that liquid, i.e.: 



We define the paired correlation factor or the radial distribution function 
g(r) by the relation: 



[ 1 . 6 ] 


P 


As we can see, this function is the ratio of the mean value of the local 
density of molecules (mean calculated at the positions, at a given time and 
over a period of time) to the volumetric density of molecules. The 
correlation factor g(r) is proportional to the probability of finding a molecule 
at a distance r + dr from another molecule. Thus, we can write the relation: 



[1.7] 


where S i . is the Kronecker delta, such that: 
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Thi s ratio [1.7] quantifies the local structure - in other words, the way in 
which the molecules are arranged in relation to one another. 


1.2.3 The curve representative of the radial distribution function 

By combining relations [1.4] and [1.6], we see that the mean number 
of molecules in the coronal volume between the spheres with radii r and 
r + dr will be: 

(d /V(r)) = Anpr 1 g(r) dr [1.8] 

In the solid crystal, only certain distances exist, and the representative 
curve for the function g(r) exhibits extremely slender peaks for these 
distances. 

In the case of the liquid, the curve representing the function g(r) has the 
shape shown in Figure 1.3. We obtain a first peak with a breadth A r/r of 
several %, which represents the distance between the first neighbors. The 
next peaks, which represent the second, third (etc.) neighbors, are heavily 
damped because of the disorder over a long distance. The function g(r ) tends 
toward 1 at a long distance, there is no longer order and therefore, on 
average, we always find the same number of molecules per unit volume as 
are present in the overall liquid. 

Figure 1.3 can be obtained by neutron diffraction or hard, very 
penetrating X-ray diffraction, such as those produced by synchrotron 
sources. 

In principle, the distribution g(r) is null for distances less than 0.5 A, 
because there is no chance of finding two molecules that close together, 
given that the order of magnitude of a molecule’s diameter is between 1 and 
3 A. Around values between 3 and 5 A, molecules may be found, and the 
local density is greater than the overall density. Thus, g(r ) is greater than 1 . 
Between the first series of neighbors and the second, there are few 
molecules, and the factor g(r) drops back below 1. The second maximum 
corresponds to the second neighbors, which are less precisely localized, and 
therefore have lower local densities - hence the damping effect seen here. 
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The more the peaks are damped, the more negligible the influence of the 
central atom. Thus, it is easy to understand that we can often content 
ourselves with the influence of the nearest neighbors, i.e. those which 
correspond to the first peak. 



Figure 1.3. Paired distribution function for a liquid 


We can see that, apart for a few exceptions, the local order and the 
intermolecular distances of the maxima are the same in the solid and 
the corresponding liquid. The peaks shown in the solid are very slender, but 
the first peak is situated at the same value of r. 

Diagrams such as Figure 1.3 are very useful, because they enable us to 
calculate two statistical values: 

- the mean distance of the first neighbors. This value is given by the first 
maximum point on the curve. The breadth of the peaks shows the variation 
of the distances around the mean value due to the ordering of the molecules 
and to their agitation; 

- the mean value of the number of first neighbors. To calculate this, we 
decide that the first neighbors are those which are found at distances 
between 0 and r min . This value is the abscissa of the minimum which follows 
the first maximum (see Figure 1.3) on the plot of g(r). Thus, for the number 
of first neighbors, we can write: 


z = j d N(r) = A7ip | r 2 g(r)dr [1.9] 

o o 

In a liquid, unlike with a crystal, this number may not necessarily be an 
integer. 
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1.2.4. Calculation of the macroscopic thermodynamic values 

We shall show that, like the partition function, the radial distribution 
function contains all the information pertaining to the thermodynamic 
definition of the liquid. Therefore, it can be used to calculate the 
macroscopic values such as the internal energy, the pressure, an equation of 
state or the heat capacities. 


On the basis of relation [1.8], we can write the differential of the internal 
energy due to the interactions in the form: 


Akt 2 pg{r)s{r)dr 
2 


[ 1 . 10 ] 


Hence, by integrating over the whole volume: 


U-U 


pf 


R T 


— — — f 4 7tr 2 g(r)e(r) d r 
2kj{ 


[ 1 . 11 ] 


U j denotes the internal energy of a fluid with no interaction, i.e. the 

molar internal energy of the perfect gas which, according to the theorem of 
equal distribution of energy, has the value: 


U 


pf 


3R T 

~Y~ 


[ 1 . 12 ] 


From this, we can deduce the internal energy: 


U _ 3 
KT~ 2 


P 


2k 


Anr 2 g(r)e(f)dr 


[1.13] 


As we have the expression of the internal energy, which is a characteristic 
function in variables V and S, we have all the necessary information to 
define the phase. 

To calculate the pressure, we need to have the differential of the internal 
energy in variables P and T, an expression which is of the form: 



d T- 


fdp ^ 


p+\ 


IdT) 

V 


larj 

V _ 


d U = T 


dV 


[1.14] 
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Then, we calculate the pressure on the basis of the derivative of the 
internal energy in relation to the volume, which gives us: 


P 

pRT 


= 1 - 


2^-n ip 

3k J 


j r 2 g(r) 

0 


d £{r) 
dr 


dr 


[1.15] 


By substituting the value of r found in relation [1.5], written for one 
mole, back into the above expression, then it is easy to write the equation of 
state: 


PV = «N a 


1- 


ip 

3k J 


j r 2 g(r) 
0 


d e{r) 
dr 


d r\RT 


[1.16] 


Similarly, the material derivative of the internal energy in relation to the 
temperature enables us to easily calculate the heat capacity at constant 
volume. 


We can now calculate all the other functions, particularly the 
compressibility at constant temperature Xt, which gives us: 


Xt 


1 

fdv' 

\ 1+N J, 

V 

U p. 



\g(r)-\]Ajrr 2 dr 

RT 


[1.17] 


Thus, we have shown that knowing the radial distribution function 
enables us to completely define the phase in thermodynamic terms. 

We know that a second way of calculating the macroscopic values is to 
use the canonical partition function. This is the method that we shall use 
from hereon in. To do so, we must construct a structure of the liquid, in 
order to be able evaluate the terms of interaction in the canonical partition 
function. Various techniques are used. We shall describe four such 
techniques: Guggenheim’s and Mie’s models, extrapolated respectively from 
the gas and solid models, the Lennard-John and Devonshire cellular model 
and the cell/vacancy model. 


1.3. Models extrapolated from gases or solids 

In light of the proximity of the structure of liquid, firstly to that of a gas 
(in terms of the mobility of the molecules and the disorder at long distance) 
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and secondly to that of a solid (in terms of the presence of an order over a 
short distance and figures of X-ray diffraction), the earliest models 
developed were extrapolations either from a model of a gas or from one of a 
solid. These models can be used to calculate the radial distribution function 
and the canonical partition function. As we know that only one of the two is 
necessary, in our discussion below, we shall restrict ourselves to calculating 
the canonical partition function. 


1 . 3 . 1 . Guggenheim’s smoothed potential model 

This model [GUG 32] is extrapolated from the imperfect gas model, 
which can be used to calculate the second coefficient of the virial (see 
section A.3.4 in Appendix 3). The canonical partition function then takes the 
form of equation [A.3.41]. 


From this, we deduce the configuration integral due to the interactions 
and to the volume, in this case the volume of slightly imperfect gases: 


I,= 


yN 

xT 


1 - 


N^T) 

V 


\ N 


[1.18] 


Using the notation v m to represent the volume per molecule (V/N), or the 
molecular volume (which must not be confused with the volume of a 
molecule), and using Stirling’s approximation [A.3.1], this expression takes 
the following equivalent form: 

/ / =exp(A)(v,„-5 AA (^)) A, [1.19] 

According to relation [A.3.40], the term B A A is a function only of the 
temperature. 

We can use such an expression for a highly-imperfect gas or a liquid, 
supposing that the term Zf A \ is also a function of the volume. The difference 
v m -B AA (T,v m ) will therefore represent the free volume per molecule v f and the 
above relation will then be written: 

I, =exp(X)(v / (vj)' V 


[ 1 . 20 ] 
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Supposing the potential function to be more or less constant, an initial 
model of the liquid state was proposed. 

Thus, we suppose that each molecule moves through a uniform field of 
potential, the lattice energy, (-£) (attractive) which will be determined by 
the mean number of near neighbors at a given distance from the molecule, 
and will essentially be a function of the number of molecules per unit 
volume, i.e. a function of the volume per molecule v m . The contribution of 
the interactions to the canonical function can thus be written by 
supplementing it with the exponential term corresponding to that uniform 
potential. Thus, we obtain: 


I i = exp(A)(v / (vj)" 


exp 


V J 


[ 1 . 21 ] 


This can also be written as: 


I,= 


exp(l)v,(vjexp 


fOO 

k B r 


[ 1 . 22 ] 


Hence, in light of relation [A.3.42], and with Stirling’s approximation 
applied, the canonical partition function for the fluid will be written thus: 


lnZ c = 7Vg(V, "' ) + Aln 
k B T 

+N + N\az^-N\aN 


(2nmk B T) V2 v f (v m ) 


V 


[1.23] 


On the basis of relations [1.23] and [A.3.48], we can calculate the 
Helmholtz energy F: 


^ = -£<rJ- k B nn 


(2^-wk B 7 T ) 3/2 v / (v m )^ 


-k B r-k B rinz int 


[1.24] 
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From this, we deduce the pressure: 


P = - 


\dV j T n 


MV m ) + k ^ V A V m) 


3v.„ 


3v„, 


[1.25] 


Most of the time, the properties of the liquid are insensitive to variations 
in pressure and it is therefore correct to consider the pressure to be zero. 
Thus, we shall have: 




3v 


9v,„ 


[1.26] 


Flence, the compressibility factor Z = PV/RT is essentially null regardless 
of the volume. 


The molecular Gibbs energy is: 

N N N 

(27rmk n T) i 2 v f 

= -e - k B T In — — k B T - k B T In z mt + Pv m 


[1.27] 


The molecular internal energy is given by: 


N 


^(F/T)\ = _ e+ 3_ 

2 B N 


V dT , 

Thus, the molecular enthalpy is: 
H U V 




— = — + P— = -e + -k n T + ^- + Pv 
N N N 2 N 

The molecular entropy is given by: 

. F , ... 


*_ = N_ 

N dT 


= k D In 


(2nmk B T) v f v m ) 


V 


5 9hnv 

+ — k R +k R r - + k R lnz. 

2 dT 


[1.28] 


[1.29] 


[1.30] 



Pure Liquids 13 


In view of the variation of the free volume with temperature 
(see section 1.3.3) and of the independence of the volume from pressure 
(Zt =0), we find: 


— = k B In 

N B 


r (2/nnkjf' 2 v f )^ 


-k B ln — + ^k B +k B lnz int [1.31] 

v... 2 


In relations [1.27], [1.28] and [1.29], the pressure can be taken to be null, 
in keeping with relation [1.26]. 


1.3.2. Mie’s harmonic oscillator model 


This time referring to the local order in a liquid similar to that of a solid, 
the potential function is given a form very similar to that of a harmonic 
oscillator. Thus, a second model of the liquid state [M1E 03] was put 
forward. Beginning with the quasi-crystalline model of a liquid, we suppose 
that each molecule is in a field of potential whose minimum is £o(v m ), and 
that the molecule moves through that field corresponding to a three- 
dimensional harmonic oscillator of frequency V, which is also a function of 
the volume per molecule v m . We use the symbol r to denote the distance 
from the center of the molecule to the center of the cavity where the 
minimum potential is in force. At that distance, the molecule would have a 

fji ( 2ttv) 2 r 

potential energy ~(e 0 + k B T) + so, if we integrate for all 

possible positions of the molecule, the configuration integral for the partition 
function is found to be: 


1 , = exp 


N 


(go + 

k J 


1 4 nr 2 exp 


^ m{l7tv) 2 r 2 ^ 

2k B r 


dr 


[1.32] 


After integration, this gives us: 


1 1 =s exp 


(go + k B 

k B r 


( k J ) 

3/2 | 

ylKmv 2 J 

J 


[1.33] 
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Using approximation [A.3.31] for the vibration, the complete canonical 
partition function is then written: 


Z c = 




Z- 


V h 2 , 


V hv j 


[1.34] 


Hence, if we take relation [1.33] into account: 


Z c = 



r In tnk B T 5 

3/2 

N | 

r r 

Z in, 

l h 2 J 


| 

^exp 


( £ o + k B Z) 

k,T 1 

3/2 'j 

k B r 

y 2nmv 2 , 

J 


[1.35] 


By switching to logarithms, we obtain: 


lnZ c = jVg ° (V " , ^ +3Aln 
k B r 


V hv/ (v m )y 


+ N + N In z. 


[1.36] 


This is the partition function of the liquid, given by Mie’s 3D harmonic 
oscillator model. 

Based on relations [1.36] and [A.3.48], we find that the Helmholtz energy 
per molecule is: 


F 

— = —£ a -3k„rin 
N ° B 


v hu J 


-k B r-k B rinz mt 


[1.37] 


In the same way as we did above, we deduce the expressions for the 
different functions: 


p = _SF = 6^_ }k y_ tav 

dV dv dv m 


= 0 


[1.38] 


Ty = ~ £ on ~ ?>K T In - k B Z - k B Tln z int + Pv m 
N hv 


[1.39] 


- = -£,+ 3k R r + -^ 
N ° B N 


[1.40] 
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— = -£•„+ 3k T + ^L + Pv m 
N ° B N 


In all these relations, the term Pv m is negligible. 


The entropy per molecule is given by: 


S_ 

N 


(-1 

-AA- = 3k B In 
dT B 


kj_ 

hv 


+ 4k B + k B In z in 


[1.41] 


[1.42] 


This expression is independent of the volume. 

Thus, we obtain two different series of expressions. We usually use the 
smoothed potential model we link the properties of gases to those of liquids, 
and the harmonic oscillator model to link the properties of liquids to those of 
solids. 


NOTE 1.1.- Relations [1.24] and [1.37] may be identical for a certain 
temperature and a certain volume per molecule, identifying e with £ 0 and 
attributing the following value for the molecular volume: 


V 


2 xmv 2 


[1.43] 


1.3.3. Determination of the free volume on the basis of the dilation and the 
compressibility 

The free volume of the liquid, which we need to know in order to exploit 
Guggenheim’s model, can be determined by a variety of methods: velocity 
of propagation of sound, vapor pressure, measurements of dilatation and 
compressibility. We have chosen to discuss this latter method. 


In view of their definitions, the volumetric dilation coefficient and 
compressibility coefficient enable us to write: 


A 

Zr 






V,N 


[1.44] 
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In light of relation [1.26] (with £ and V/ being functions only of v m ), we 
can write: 


' <)P \ 3 In Vj 

^ t )v,n B 


[1.45] 


However, if we consider that the molecules are arranged in a cubic lattice 
with centered faces, spaced the length a apart, we can take the following 
value for the free volume for the molecules: 


4 K 

T 


( a-Df 


[1.46] 


Using relation [1.51], which we shall demonstrate later on (see relations 
[1.50] and [1.51] in section 1.4), we can write: 


3 In v f 
3v 

m 



[1.47] 


By substituting this value back into equations [1.44] and [1.45], we find 
the value of the molecular free volume: 


V/ = 


Wgkj 'Izl 

3v;,/? 3 


[1.48] 


If, for P and % T , we take orders of magnitude of, respectively, 10‘ 3 K _1 and 
10' 4 atm _1 , we obtain, e.g. for chloroform: 

N a v/=0.44 cnrVmole [1-49] 

Thi s value is perfectly acceptable. 


1.4. Lennard-Jones and Devonshire cellular model 

This model [LEN 37] is based on Figure 1.4. Each molecule is inside a 
spherical cage - the cell - whose radius is a. This sphere is the molecule’s 
mean sphere of influence. 
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The energy of interaction between two molecules is expressed by 
Lennard-Jones’ relation [A.3.44]. This interaction is limited to the 
molecule’s z nearest neighbors. 



This number z is the coordination index linked to the cell geometry. We 
shall suppose that the molecules occupy the sites of a cubic lattice with 
centered faces, and therefore the coordination index is z = 12. 

The volume of liquid is divided into cells centered on each molecule, 
whose near neighbors occupy the medium from the vertices of a cube with 
side length 2 a / -Jl . Each molecule which is a near neighbor of the original 
one thus belongs to four cells, and each cell contains 1 + 12/4 = 4 molecules. 
Hence, the volume of the cell is such that: 

4v m =8n 3 /(V2) 3 [1.50] 

and therefore: 

« 3 =v m V 2 [1.51] 


The translational canonical partition function with interactions can be 
written, on the basis of expression [A.3.38], taking account only of the z 
molecules that are near neighbors of each molecule i. 


^ C(<) N\ 
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[1.52] 
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On the basis of relation [A. 3. 24], using e (co N ) to denote the double sum 


appearing in the previous relation, the canonical partition function becomes: 




Figure 1.5. Distance between a molecule and one of its near neighbors 


We shall now evaluate the double sum e (of ) . The near neighbors of a 

molecule i are, in fact, situated at different distances dj from it, and those 
distances change over time. Rather than moving both the molecule i and its 
near neighbors at once, we shall suppose that the center of the cell is 
stationary, that the molecule i moves in concentric circles of radius r around 
that center, and that the near neighbors are affixed to a concentric sphere 
with radius a. The radius r varies between 0 and a. Consider the plane 
passing through the molecule i at point I, one of its near neighbors at M and 
the center of the cell O (Figure 1.5). 

In order to simplify the multiple integration appearing in expression 
[1.53], we shall create spherical symmetry and average the distance from the 
atom i along a radius r to all its near neighbors. The energy £ (o) N ) can be 


written: 



[1.54] 
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( £ (d)) is the mean mutual energy between the atom i and a near 

neighbor when the atom i is on the disc of radius r. This energy is a function 
of the distance, which is given by the following equation (see Figure 1.5): 

d = -Jr 2 + a 2 + 2 ar cos 0 [1.55] 

The mean energy is: 

K 

(e{d)) = ^£{d)sm0&6 [1.56] 

o 


Using relation [A.3.44] for the potential energy of interaction and 
substituting into it the value given by relation [1.55], we find: 
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with the following meaning for the functions “1” and “m” that appear in 
equation [1.57]: 


i 


f 


1 + 12^ + 50 


( 1 it 


2 A 


+12 


a 

d 2 


\ a j 


A J2 V f d 2 V 

+ ' 

V a j j 


i-Y 


-l 


[1.58] 



' d 

1 + — 

v a j 


V 


d lxA 
" J 


-1 


[1.59] 


Figure 1.6 shows two curves illustrating the variations of our potential 
energy as a function of the ratio d/a for two values of the ratio dja equal to 
0.942 (part a) and 0.681 (part b). A posteriori, these two functions provide a 
justification: the first, for the approximation of the smoothed potential 
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theory, as the potential energy is practically constant (part a), and the 
second, for the approximation of the 3D oscillator theory, with the potential 
energy showing a near-parabolic shape as a function of the distance (part b). 


► 



Figure 1.6. Potential for interaction of a molecule in a liquid according to 
Lennard-Jones and Devonshire, a) for dfa = 0.942; b) dfa = 0.681 


By substituting expression [1.57] back into relation [1.53], the 
contribution of the translational motion to the canonical partition function 
can be written as: 


^ c(0 ]\f\ 
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[1.60] 


According to Lennard-Jones and Devonshire, the integration limit in 
equation [1.60] is of no importance, because the greater part of the 
contribution is made by small distances, particularly where d < all. 


If we set x = d/a, then the logarithm of the translational partition function 


is: 
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We have chosen z = 12, and tj(£ 0 ,d 0 ) denotes the function defined by: 


n(£ 0 ,r 0 )= j * 1/2 exp 


k B r 


f d V2 
u 0 


V a r) 

-2 


f d ^ 

u o 


v a rj 


10 ) 

m(x) 


dx 


[1.62] 


c/q/ Cl 

^ = -9 
k J 

-5l = -10 

k J 

0.942 

0.00180 

0.00161 

0.918 

0.00295 

0.00269 

0.891 

0.00515 

0.00478 

0.858 

0.00964 

0.00916 

0.818 

0.01957 

0.01920 

0.765 

0.0437 

0.0445 

0.730 

0.0676 

0.0700 

0.681 

0.1069 

0.1125 


Table 1.1. Values of the function i)(e 0 ,d 0 ) 


The energy £{a) is given by: 


e(a) = e 0 



+ 



12 


V “ ) V a ) 


[1.63] 


Table 1.1 gives a few values, which are easy to calculate automatically, 
for this function for two values of the ratio fo/'k rT and different values of the 
ratio dj a. 

Thus, if we accept the hypothesis of a cubic stack with centered faces, 
i.e. y=\f2 , and if we know the molecular volume, the translational partition 
function contains only two parameters linked to the substance: d 0 and 
which are the two parameters that play a part in the expression [A. 3. 44] (in 
Appendix 3) of the Lennard-Jones interaction potential. 

In order to compare the result to other models and to experimental results, 
we need to deduce the expressions of the thermodynamic functions. From 
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relation [A.3.48], we deduce the expression of the molecular Helmholtz 
energy function: 
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-k B r In ( In a 1, ) - k B T In z int - k B T In 2y/2n7 ](£ 0 , d 0 )v m 


Here, z int represents the contribution of all other internal motions of the 
molecule to the molecular partition function (rotations, vibrations, electronic 
and nuclear spin motions). For atomic liquids, this term can be taken as 
being equal to 1 . 


Based on the Helmholtz energy, it is easy to obtain the other 
thermodynamic functions such as: 
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Thus: 
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7j,(£ 0 ,d 0 ) and T) m (e 0 , d 0 ) being two functions, such as Tj(£ 0 ,d 0 ) , of the 
two variables do and £o. They can be calculated numerically using the 
relations: 
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nS e o’ r o)= j *' 2 m(x)exp 


1/2 


12£ 0 V “/ / 



dx 


[ 1 . 68 ] 


o 
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dh>/ 42 k B T 



Figure 1.7. Isotherms calculated using the Lennard-Jones and Devonshire model 


NOTE 1.2.— In expression [1.60], by comparison with the translational 
canonical partition function for a perfect gas (relation [A.3.26]), we can 
define the free volume of the molecules as: 



[1.69] 


= 2mtie Q ,d 0 )Yv m = 2 ^ 2 rni(e 0 ,d 0 )v m 

Figure 1.7 shows a few forms of isotherms in the representation 
dlP/JlkJ as a function of sf2v m / d-l . Curve (III) is that of a perfect gas, 
curve (II) is obtained for 12£ 0 / k B T = -9 , and curve (I) for 
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1 2e 0 / k B T = -10. It is noteworthy that curve (II) seems very similar to the 
critical isotherm which is given by: 

T c=~ [L70] 

3k B 

The values obtained for certain liquids (see Table 1.2) show a satisfactory 
degree of accordance with their experimental values. 


Substances 

d 0 ( A) 

—£ 0 (10‘ 22 joules/molecule) 

T c calculated 

T c observed 

H 2 

35.3 

4.25 

41 

33 

Ne 

29.2 

4.89 

47 

44.47 

n 2 

72.5 

13.25 

128 

126 

A 

56.2 

16.5 

160 

150.66 


Table 1.2. Values of the critical temperature, found experimentally and calculated 
by the Lennard-Jones and Devonshire model 


It is a fairly laborious task to rigorously calculate the critical volume, but 
from Figure 1.7, it seems we can choose the critical volume such that: 


= 2 so that v =dls ! 2 
d 

u o 


[1.71] 


We can see that this value is far too low. Indeed, it yields a value of 0.7 
for d^P / \fik B T , instead of 0.3, which is the result found experimentally. 

Thus, the Lennard-Jones and Devonshire cellular model can be used to 
calculate thermodynamic functions with only two adjustable parameters. In 
section 1.7, however, we shall demonstrate that the results obtained are very 
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approximate, so there is a necessity to perfect the model. This was the 
purpose of the cellular and vacancies model developed by Ono and Eyring. 


1.5. Cellular and vacancies model 

Whilst it does represent real progress in relation to the previous two 
models, the Lennard-Jones and Devonshire model discussed above has a 
serious shortcoming-it is incapable of taking account of two dynamic 
properties of liquids: the phenomena of viscosity and self-diffusion. In order 
to take account of these properties, Ono [ONO 47] introduced the concept of 
vacancies, comparable to that which takes account of conductivity and 
diffusion in the solid phase. Ono considers that certain sites in the pseudo- 
lattice, or if you prefer, certain cells described in the above model, are not 
occupied, forming what we call vacancies. Thus, on average, over time, a 
molecule i will be surrounded by z, first neighbors in accordance with: 

z i = y i z [1.72] 


y>i appears as the fraction of first-neighbor sites occupied around the 
molecule i. Therefore, v, is a short-distance order index, whose value is zero 
when the central molecule i is surrounded only by holes (i.e. no molecules), 
and 1 if all the cells neighboring the central molecule are occupied (see 
section 3.2.1). Its spatial mean would be <y,>, and would correspond to the 
mean value of the number of first neighbors <z,> determined by the first 
maximum of the radial distribution function demonstrated by X-ray 
diffraction. The number z, which is the coordination index, is in fact that 
maximum possible number of first neighbors, given by the chosen structure; 
that value is often taken to be equal to 12 for cubic cells with centered faces. 

In order to take account, individually, of the environs of each molecule, 
we divide the liquid volume V into L cells ( L>N ) with respective volumes 
T\, Ti ... T\ . The configuration integral, which is expressed over the whole of 
the volume V, and plays a part in relation [1.53], will be replaced by a sum 
of partial integrals, each of which corresponds to an individual cell. Relation 
[1.53] then becomes: 
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The sum contains L N terms which correspond to L different cells, in 
relation to which the coordinates of the different molecules are expressed. 
Because each cell is supposed to be sufficiently small to contain at most one 
molecule, and sufficiently large so that the intermolecular forces are 
practically no longer felt beyond the immediately adjacent cells, the sums in 
expression [1.73] contain only L\I(L-N)\ non-null terms. 

The energy e[co N ^ in expression [1.71] will therefore be rewritten, 
instead of expression [1.54], in the form: 

*K) = Et e{a) + ^\_{£{d i ,y i ))-£{a)] [1.74] 

!=1 ^ (=1 ^ 


The distance d l is given, for each cell, by a relation similar to 
expression [1.55], namely: 

d t = yji] 2 + a 2 + 2ar f cos 0 [1-75] 


We use the notation v c to denote the volume of a given cell, and for the 
cubic lattice with centered face, we have: 



[1.76] 


Note that this cellular volume differs from the molecular volume 
v m = V/N, because we no longer always have a molecule in each cell. 

NOTE 1.3.- The volumes v c and v m have a known ratio, because we have: 


v,„ 



[1.77] 


Hence, instead of relation [1.61], the translational canonical partition 
function is written as: 


Z = 

r(t\ — 
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2k B r 
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[1.78] 
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The function j (y ,) is homogeneous with a volume. It is defined by: 



[1.79] 


If all the sites are occupied, v, = 1 and j(y,) is identified with the free 
volume logarithm v/ from relation [1.69]. Ify, = 0, then j(0) is the logarithm 
of the volume of the cell v c . 

We can see that the function j(y,) is not a simple expression, and the 
various authors have been led to simplify it by a linear form of y, as a 
function of the logarithm of Vj and v c so as to satisfy the boundary values of 
j(y j). Thus, Ono proposed the expression: 



[1.80] 



Figure 1.8. Potential energy curve for a molecule occupying a 
more favorable position than a neighboring vacancy (from [REE 64]) 


Eyring and his collaborators [REE 64] put forward the formula: 


j O, ) = y, In ( v f g, ) + ( 1 - y, ■ ) In v c 


[1.81] 
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Eyring’s function differs from Ono’s only by the introduction of a factor 
g, known as the degeneration factor, which is introduced because it is 
natural to suppose that, of all the places that are available, on average, for a 
molecule, some are more favorable than others, energetically speaking, 
simply because of the organization of the other molecules. This difference in 
energy between the most probable place and a less-probable near neighbor 
(see Figure 1.8) must be proportional to the interaction energy zs(a) / 2 , and 

inversely proportional to the number of vacancies n h = z(l- (>’,)) ■ Thus, 
this energy difference would be of the form: 


As — 


kzs(a) 
2 n h 


ks(a ) 

T-W) 


[1.82] 


k is an adjustable constant of proportionality. Thus, the degeneration 
factor due to the vacancies present around a molecule would be written as: 


g = l + z(l-(y,.))exp 


As 

k J 
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[1.83] 


The authors show that if we take account of relation [1.81], using the 
Bragg-Williams approximation (A e= 0 in g, see section 3.1.2) and Stirling’s 
approximation [A.3.1], the translational canonical partition function [1.78] 
assumes the form: 
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Using this expression, the thermodynamic functions - particularly the 
Helmholtz energy - can be calculated. 


1.6. Eyring’s semi-microscopic formulation of the vacancy model 

The expressions used by Eyring in the vacancy model become 
complicated and tricky to calculate numerically. In order to remedy this 
situation, Eyring and his collaborators [EYR 6 1 ] adapted the vacancy model 
to a semi-microscopic model, by replacing Lennard-Jones’ interaction 
functions with macroscopic values, involving the molar volume of the solid 
v Wsoi) at its melting point and that of the liquid \> Mliq) at the temperature of 
study. Observing (except for a few very rare cases, one of which is water, 
which is not a non-associated liquid) a significant increase (around a twofold 
increase) in the molar volume upon transitioning from the solid state to the 
liquid state, the authors model a liquid as being a two-component solution: 

- molecules, which behave like a molecule in a solid, i.e. three 
vibrational degrees of freedom. This is the model of the short-distance lattice 
aspect; 

- vacancies, which behave like a gas, and therefore have three 
translational degrees of freedom, which we shall suppose to be perfect, with 
non- localized objects that are free to move around, which will create 
disorder over a long distance and mobility of the species. 

Of course, the movement of a vacancy is, in fact, the movement of a 
molecule neighboring that vacancy. 

Solids, just below their melting point, are assumed to contain no 
vacancies. If there are any, they are few in number in relation to the 
molecules. In a liquid, on the other hand, the number of vacancies is of 
the same order of magnitude as the number of molecules. If N denotes the 
number of molecules behaving like a solid, the number of vacancies would 
be: 



30 Modeling of Liquid Phases 


Thus, the total number of cells would be: 


N c =N l +N = N 


v °(«?) _ v 0 (sol) 


[ 1 . 86 ] 


Hence, the fraction of sites with molecules would be: 
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and the fraction of sites with vacancies would be: 
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[ 1 . 88 ] 


The mean number of vacancies neighboring a molecule (z is the 
coordination index of the lattice) would be: 


(«i) = z — = ■ 


0 (sol) 


Mliq) 


[1.89] 


By applying relation [A.3.36] to both components of the solution, we can 
calculate the canonical partition function on the basis of that of the localized 
molecules and non-localized vacancies, so that: 
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As the vacancies behave like a perfect gas with three translational degrees 
of freedom, we have: 
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For molecules which behave like a solid, with three vibrational degrees of 
freedom, if we ignore the residual vibration, we have: 
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g 


1-exp 


' hv^ 

V J 


[1.92] 


The degeneration coefficient is calculated in a similar manner to that used 
to obtain relation [1.85]. Its value is: 


g = 1 + n exp 


As 

k7 


[1.93] 


If A S U is the Helmholtz energy of sublimation of the solid, the variation 
in energy As (Figure 1.8) will be proportional to the sublimation energy and 
inversely proportional to the number of vacancies Nl- We can write this in a 
similar manner to expression [1.82]: 
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The partition function of the molecules with solid behavior would 
therefore be: 
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Using expressions [1.90], [1.91] and [1.95] for the overall partition 
function, and after application of Stirling’s approximation, we find: 


Nv 0i 


sol ) 

lnZ r =- n ... . 

L ,,0(hq) 


AU 


R T 


-3 In 


1-exp 


f hv ^ 


V 


+ 


In 


l + z 


yOm _ yO(SOl) 


v 


N(v° m -V 0(so,) ) 


0 (liq) 


„0 (//?) 


— 31n- 


-exp- 


kAUv 0(so,) 


RT{ 


y 0 (liq) _ yO(so/) ' 

0(/i?) 




[1.96] 


(2/rmk B T) 


- + ln 


ev 


N 


We can see that the canonical partition function, and therefore all the 
thermodynamic functions (particularly the Helmholtz energy) depend only 
on the single adjustable parameter k defined in expression [1.94]. Hence, this 
model is a simple and powerful tool. 

The model we have just looked at is that which applies to atomic liquids, 
such as argon, for instance. Eyring and his collaborators carried out parallel 
tests, applied to the cases of molten salts and liquid metals. 

NOTE 1.4 - It is worth noting that Eyring’s model, for the partition function 
(relation [1.90]) is tantamount to mixing Guggenheim’s (section 1.3.1) and 
Mie’s (section 1.3.2) models. These models gave a good account of the 
properties of liquids respectively in the vicinity of a gaseous state and of a 
solid state. 


1.7. Comparison between the different microscopic models and 
experimental results 

A variety of comparisons have been offered by the different creators of 
models: comparisons between a model and the experimental results, 
comparisons between different models, and comparisons between results 
produced by a model and those produced by simulation calculations. Indeed, 
the calculation methods used for statistic simulation lend themselves very 
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well to this type of problem. Examples include the static Monte Carlo 
method, based on equation [1.3], or the molecular dynamics method, based 
on the fundamental law of dynamics (see Appendix 1). 

In terms of comparison with experimental values, we shall give the 
example of the variation in heat capacity at constant volume as a function of 
the temperature calculated by Eyring’s semi-microscopic method. 
Remember that it is a cellular model including vacancies and a degeneration 
coefficient (see section 1.6). Figure 1.9 illustrates such a comparison and 
exhibits good accordance between the results obtained by the model and the 
experimental results. 



ft 



4 - o Experiment 


8 ) 90 100110 120 130 14 ^ ^° K) 


Figure 1.9. Comparison between the observed values of 
the heat capacity at constant volume and those calculated using 
the cellular and vacancies model by Eyring et al. [EYR 61] 


Figure 1.10 shows the comparison of the result of the same model of the 
radial distribution function curve for argon at a temperature of 84.4K, 
against the experimental result. Once again, we see excellent accordance. 

Certain comparisons are made between the measured values and those 
calculated by a model, for the critical values - particularly the critical 
temperature and critical pressure, using the conditions: 
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Figure 1.10. Comparison between the experimental value and that calculated by Eyring ’s 
model for the radial distribution function of Argon at 84.4K ( data from [YOO 81]) 


Substance 

C(°K) 

P c (atm) 

Calculated 

Observed 

Calculated 

Observed 

Neon 

55.41 

44.47 

37.65 

26.86 

Argon 

154.44 

150.66 

58.72 

48.00 

Krypton 

208.33 

210.60 

69.68 

54.24 

Xenon 

287.80 

289.80 

74.89 

58.20 


Table 1.3. Comparison of the observed values and those calculated by the Eyring 
model, for the temperature and the critical pressure (data from [EYR 58]) 


Certain data appear in Table 1.2 for Lennard-Jones and Devonshire’s 
model (see section 1.4). Others are given for the solids of rare gases in 
Table 1.3, and pertain to Eyring’s model (see section 1.5). 

Note that both models yield satisfactory results on this point. However, it 
is important to apply the comparison to several types of results. For example, 
Figure 1.11 shows that, for the representation of the distribution function, 
Fennard-Jones and Devonshire’s model, Eyring’s model and the calculations 
performed by numerical simulation are very similar. Meanwhile, 
Figure 1.12, which gives the variation of the compressibility coefficient as a 
function of a reduced volume, illustrates the significant behavioral difference 
between the molecular dynamics simulation and Eyring’s model, on the one 
hand, and Fennard-Jones/Devonshire’s, Guggenheim’s (see section 1.3.1) 
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and Mie’s model (see section 1.3.2), on the other. It seems that the important 
point which divides Eyring’s model from that of Lennard- Jones and 
Devonshire is more the introduction of the degeneration coefficient than the 
variability of the number of molecules that are near neighbors of a given 
molecule (z ; ). 



Figure 1.11. Comparison of the curve of the radial distribution function between the 
calculations of molecular dynamics and different models (according to [YOO 81]) 


PV/RT 



Figure 1.12. Comparison of the results obtained on 
the compressibility factor ( data from [REE 64]) 
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Thus, we can see that it is important to examine the validity of a model 
by comparing several results produced by that model. Additionally, a good 
model of the structure of liquids must also satisfy the interpretation of 
properties other than the mere thermodynamic values that are of interest to 
us here, e.g. surface tension, viscosity and self-diffusion. The major 
advantage of Eyring’s cellular and vacancy model with a degeneration 
coefficient is that it also takes account of the dynamic properties of liquids. 



2 


Macroscopic Modeling of 
Liquid Molecular Solutions 


The representation of solutions must give the expressions either for the 
molar excess Gibbs energy or the activity coefficients, using a given 
convention for each component as a function of the variables of the problem, 
which are usually the pressure, temperature and the composition of the 
solution. In this chapter, we shall only discuss the case of liquid molecular 
solutions. In such solutions, usually, the pressure is no longer a variable of 
the system. 

A distinction must be drawn between macroscopic representations and 
simulations using a microscopic model. 

The advantage of a representation is that we have an expression for the 
excess Gibbs energy or the activity coefficients, as a function of the 
solution’s composition and possibly of the temperature, based on experience, 
which can be used in expressing the law of action of the masses or in 
constructing phase diagrams. 

The advantage of a microscopic model is that it shows the same qualities 
as the representation, with an added predictive element and a physical 
meaning for the different properties that are involved. 

We often work with binary solutions, which are obviously simpler, but 
significant efforts are being invested in modeling solutions with more than 
two components, based on the knowledge of binary solutions formed with 
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different pairs of components of the solution. Another interesting method is 
the modeling of binary solutions using the functional groups of the 
molecules. 

In this chapter, we shall discuss macroscopic models. The next chapter 
(Chapter 3) will describe microscopic models. 

The expressions of the activity coefficients, given to represent a solution, 
must be self-consistent, and in particular, they must satisfy relation [A.2.23], 
which is derived from the Gibbs-Duhem relation. 


2.1. Macroscopic modeling of the Margules expansion 

A helpful and simple representation of a non-perfect solution is given by 
a limited expansion. 

Consider a solution with two components, 1 and 2. It is always possible 
to represent it by expanding the logarithm of the activity coefficient for 
component 1 in reference (I), in a MacLaurin series of the molar fraction x 2 
of the other component - for example: 

T In y[ * — .1 + (X .v. + .w + S i x 1 + £ l x 2 +... [2.1] 

and, absolutely symmetrically for the second component: 

Tin Y~! 1 = A 2 + (X 2 x i + /? 9 x 2 + 8~,x 2 + £~,x^ + ... [2.2] 

These expansions are known as Margules expansions. 

Conventionally, the activity coefficient of a component tends toward 1 if 
the molar fraction of the other component tends toward zero, so the limited 
expansion does not contain a constant term. 

4=4=0 [2.3] 

In addition, according to relation [A.2.23], we should have: 


jCj d In y l + x 2 d In y 2 = 0 


[2.4] 
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From this, we can easily deduce the following relations, if we examine 
the first five terms in the expansion: 


a x - a 2 - 0 


[ 2 . 5 ] 

A = A + fi + A and A 

li 

|3S> 
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I- 05 
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[ 2 . 6 ] 
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-T = — 1-2-2- and -1 = 

L-2-L 

4 4 3 4 

4 3 

3=4 


[ 2 . 8 ] 


Thus, the Margules expansion contains no constants or first-degree terms. 
The first non-zero term is the second-degree one. 

It is worth noting that, if we look at an expansion with two non-null 
coefficients [i. and 8 j , we can write, for component 1, for instance: 

TiU fi = /3 l +S ] x 2 [ 2 . 9 ] 

x 2 

We can attempt to plot a straight line on the experimental points in a 
representation T\nyJ x' 2 as a function of x 2 . If the representation is 
properly linear, we can easily determine the two coefficients. The ordinate at 
the origin of this straight line is and its slope is d x . 


2.2. General representation of a solution with several components 

Still with the aim of having mathematical expressions for the 
representation of the solution, Redlich and Kister offered a representation 
that provides an expansion of the excess Gibbs energy, a pure-substance 
reference in the same state of segregation as the solution (reference (I)), the 
equivalent of the Margules expansion for the activity coefficients. For a 
two-component solution, the polynomial expansion up to order m is written: 

m 

G:=x 1 x 2 Y j L«\x ,- xj' 

k'=0 


[ 2 . 10 ] 
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The coefficients are parameters which reflect atomic interactions. 
Because the exponents k ’ are odd and even integers, it is internationally 
agreed to write the differences (xj — x 2 ) in alphabetical order of the chemical 
symbols of the components, written as subscripts next to the molar fractions. 

For a solution containing a given number of components, relation [2.10] 
becomes: 


g:: =T, x > x iY, L T( x i- x j) k 


•*j 


The alphabetical convention applies to all the couples i,j. 


[ 2 . 11 ] 


2.3. Macroscopic modeling of the Wagner expansions 

The Wagner expansion is another form of limited expansion, this time 
used to represent the logarithms of the activity coefficients for the solutes in 
a solvent, as reference (II) - an infmitely-dilute solution. 


2.3.1. Definition of the Wagner interaction coefficients 

Consider a solution containing the components 1, 2, ..., i, ...j, .... 
Suppose that this solution is diluted in solvent 1. In order to express the 
variations of the activity coefficients for the solutes as a function of 
the molar fractions x, (i f 1 ), which are much smaller than 1 , let us take the 
expansion of In yf into a Taylor series in the vicinity of 1 : 


In r" =ln 7i +Y, 


^in r! 1 ^ 


j = 2 


dx t 


+ . 


J J 


[ 2 . 12 ] 


If, in that expansion, we content ourselves with the first two terms, we 
can write: 


ln^' =ln^" H + X-V/ 

J = 2 


[2.13] 
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is the value of the activity coefficient of the solute i if the solution 
is very dilute - i.e. when all the molar fractions Xj of the other solutes are 
null. 


£ / is called the Wagner coefficient of interaction of solute i with solute j. 
Thus, it is defined as: 


£ = 


din f 1 ^ 

V dx J / 


[2.14] 


By simply applying relation [A.2.27] (see Appendix 2), we can show that 
the matrix of Wagner coefficients is symmetrical, so: 

£/=£) [2-15] 

We can show, by application of relation [A.2.23] (see Appendix 2), that 
the activity coefficient of the solvent obeys the relation: 

[2.16] 


This Wagnerian representation is widely used for metal alloys with low 
degrees of addition of the alloyed elements. 


2.3.2. Example of a ternary solution: experimental determination of 
Wagner’s interaction coefficients 

Consider the ternary system constituted by solvent 1 and the two solutes 
2 and 3. According to relation [2.16], we have: 


In = In ^ °° + x 2 e\ + x 3 £ 2 

[2.17] 

In yf = In y : \ n )o ° + x 2 e\ + x 3 e 3 

[2.18] 


We can distinguish two types of Wagner coefficients: 
- the diagonal coefficients el and £ 3 ; 
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- the symmetrical coefficients, which are equal according to 
expression [2.15]: e\ =e] . 


Let us first determine the diagonal terms - e.g. e \ . By definition: 


0 In 

dx 1 


a\ 


[2.19] 


As this term does not depend on the presence of component 3, we can 
operate in the binary system 1-2, and simply measure the tangent to the 
origin of the curve giving In y" as a function of x 2 in that binary system. The 
slope of the tangent is s\ . 


The same process can be used to obtain the other symmetrical coefficient, 
this time considering the binary system 1-3, so £\ . 



Figure 2.1. Obtaining of the Wagner interaction 
coefficients for a tertiary system 


Let us now examine the case of the mixed term e \ . By definition, we 
have: 


4 = 




01n y^ 

V 3x 3 j 


[ 2 . 20 ] 


For different values of x$, we plot (Figure 2.1(a)) the experimental curve 
of the variations of In y" as a function of X 2 . By extrapolation to %2 = 0, we 
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obtain the values of In y" at x 3 — 0 for different values of X 3 . We report 
(Figure 2.1(b)) those values as a function of x 3 , and by extrapolation, we 
obtain the tangent to the origin - i.e. e\ . 


2.4. Dilute ideal solutions 

Below, we shall examine a category of solution which is a limit category, 
known as dilute ideal solutions. These solutions must not be confused with 
perfect solutions, although they do have a certain number of properties in 
common. For this reason, we shall never use the term “ideal solution” to 
speak of a perfect solution, as some authors do. 


2.4.1. Thermodynamic definition of a dilute ideal solution 

If, in a solution, all the components bar one - the solvent - are very 
dilute, the number of molecules of solute is so low that they have no effect 
on the behavior of the solvent, which could be said to exhibit perfect 
behavior. In this case, we say that we have a dilute ideal solution. Because of 
the distinction between solvent and solute, we shall use convention (II) - the 
infinitely dilute solution reference - to express the activity coefficients. 

We know that for the solvent, conventions (I) and (II) are identical, and 
as, hypothetically, it should exhibit perfect behavior, we will have: 

/o 7) = ?T=l [2-21] 

Let us apply relation [A.2.23] (see Appendix 2) to the set of components: 
x 0 dln// ) +Xx s dln^ //) =0 [2.22] 

s 

From this, we deduce that y i s II) = constant . Yet we know, because of the 
convention chosen, that if the solution is in its reference state, yf H) = 1 . 
Hence, the value of the constant is 1, and the activity coefficients for all the 
components will be f" ) = 1 at any temperature. 
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2 . 4 . 2 . Activity coefficients of a component with a pure-substance reference 

With regard to the solvent, we already know its activity coefficient in 
convention (I) (pure-substance reference), through expression [2.21], 
Similarly as above, by applying relation [A.2.23] (see Appendix 2), we 
deduce that for all solutes, the activity coefficient in reference (I) is a 
constant. This constant is independent of the composition, but its value is 
not 1 because the reference state does not overlap the domain of 
definition of reference (I). However, given that, in reference (II), the 
activity coefficient for the solute is 1, we can apply relation [A.2.21] 
(see Appendix 2). From this, we deduce that the constant is Henry’s 
constant K sII : 

Y? = K sH [2.23] 

It is only in the particular case of a perfect solution that this constant has 
a value of 1 . 


Thus, the perfect solution appears to be a special case of the ideal dilute 
solution - for which Henry’s constant is equal to 1 - but the two must not be 
confused. 


2 . 4 . 3 . Excess Gibbs energy of an ideal dilute solution 

According to expression [A.2.31] (see Table A.2.2 in Appendix 2), we 
immediately obtain: 

G; = R r£x s In /P = RT^X, In K sII [2.24] 


This excess Gibbs energy is non-null, but in that excess Gibbs energy, 
there is no contribution from the solvent ( f 0 I) = 1 ). 

We can see that because the excess Gibbs energy is not null, we have 
another reason to avoid confusing an ideal dilute solution with a perfect 
solution. 



Macroscopic Modeling of Liquid Molecular Solutions 45 


2.4.4. Enthalpy of mixing for an ideal dilute solution 

In view of the definition of an ideal dilute solution, we can write that: 


[2.25] 


Thus, given relation [A.2.24] (see Appendix 2), we have: 


h.=h: 


[2.26] 


where, for the solvent: 


H=hl 


[2.27] 


The partial molar enthalpy of the solvent is identical to its molar enthalpy 
in the pure state in the same state of aggregation as the solution. 

and for the solutes: 


9 In '/P _ h° - H s _ h° - H s 


dT R T 2 


[2.28] 


The difference: 


h° - h: = A, 


[2.29] 


does not depend on the composition. The molar enthalpy of mixing 
therefore, according to equation [A.2.17] (see Appendix 2), is: 


H? =2 >,4 


[2.30] 


and the corresponding molar value is: 


H mx = A, 


[2.31] 
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Thus, the constant A s appears as the partial molar enthalpy of mixing of 
the solute s in the solution, which enables us to write the following, in light 
of relation [2.29]: 

Jr = h] - lrf x [2.32] 


2.4.5. Excess entropy of a dilute ideal solution 

Given equations [2.30] and [2.24], we can write: 



s s 

T 


[2.33] 


This function, which does not contain any terms pertaining to the solvent 
either, is non null. Here, again, it is important to avoid confusing the 
representation of the ideal dilute solution with that of the perfect solution. 


2.4.6. Molar heat capacity of an ideal dilute solution at constant pressure 

If we neglect the variations of the enthalpy values with temperature 
(which is often acceptable within a reasonable temperature range), we have: 


s~lXS 

— 


3 H mix 

3 T 


= 0 


[2.34] 


Thus, the excess heat capacity at constant pressure of dilute ideal 
solutions is practically non-existent, which means that - like perfect 
solutions - dilute ideal solutions obey Kopp’s law of additivity given by 
relation [A.2.14] in Table A.2.1 from Appendix 2. 


2.5. Associated solutions 

A number of solutions exhibit significant differences from perfection, due 
largely to the forces of attraction that are exerted between certain species in 
solution. The idea of the model of an associated solution is to average the 
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properties of the solution by replacing the initial components with classes of 
components that include the initial components but also associated species 
(dimers, trimers, etc., associations between different molecules, or 
dissociated species), resulting from the condensation of the initial species, at 
equilibrium with them and with one another. We then consider that the new 
solution thus defined has perfect behavior. 

The associated and dissociated species can be revealed experimentally, 
because the forces of attraction are sufficiently strong to be detected by 
certain techniques, such as infrared spectroscopy. 

This method is commonly used, for example, for the study of so-called 
“weak” electrolytes, where the dissociation of the molecules into ions is 
considered to be incomplete, and there is thought to be an equilibrium 
between the ions and the non-dissociated molecules. This model is also at 
the root of the quasi-chemical method of dealing with solutions (see 
Chapter 3). 

The superposition of one or more equilibrium states between the newly- 
defined species enables us to calculate the activity coefficients or the partial 
molar Gibbs energy values for the initial components. 


2 . 5 . 1 . Example of the study of an associated solution 

In order to illustrate the model of the associated solution, we shall take 
the example of a binary solution obtained by mixing n A moles of a 
component A and « B moles of another component B, and suppose that the 
forces of attraction between two molecules of A are strong in comparison to 
k B T , and in comparison to those that are exerted between two molecules of 
B, and between a molecule of A and a molecule of B. This property leads us 
to a second description of the solution, in which we accept that two 
molecules of A can come together to yield a dimer molecule A 2 , in 
equilibrium with the monomer A. 

We describe the solution under examination from two different points of 


view: 
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One option is to consider it to be a non-perfect solution, with n A moles of 
A and n B moles de B (initial quantities), with the molar fractions x A and x B 
respectively, given by: 

x A = — — — and x B = — — — = \-x A [2.35] 

n A + n A n A + n A 

The chemical potentials will be ju A and jU B , and the activity coefficients 
y A and y B . This is the first description of the solution. 

Alternatively, we can consider the solution to be a perfect solution, with 
n \ moles of monomer A, n \ 2 moles of dimer A 2 and n ' B moles of B with 

the molar fractions: 

X ' A - ”'a ;x , a _ 

n \ +n \ 2 + n \ A: n \ +n \ 2 + n \ 

and x' B = — [2.36] 

n\+n\+n\ 

The chemical potentials will be /u' A jU' A and p' B . 

In this solution, we have a state of chemical equilibrium: 

2A = A 2 [2.R1] 

This imposes the condition on the affinity: 

A = 2jU\-jU\ 2 =0 [2.37] 

Furthermore, the law of conservation of the elements imposes the 

following relations between the amounts: 

n B =n' B [2.38] 

n A =n\ + 2n\ 2 


[2.39] 
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Using the symbol a to denote the extent of the equilibrium, we shall write 
that the amounts of the species containing A in the second solution will be: 


\=n A (\-a) 

[2.40] 

a, = n A a / 2 

[2.41] 


The sum of equations [2.40] and [2.41] does indeed produce 
relation [2.39]. 

This is the second description of the solution. 


NOTE 2.1- between the two descriptions, the molar fractions of B are not 
the same in the two cases, although the numbers n H and n\ are equal. 


2 . 5 . 2 . Relations between the chemical potentials of the associated solution 

We shall now establish a relation between the chemical potentials of the 
different components of the second description as a function of those from 
the first description. 

The Gibbs energy of the solution does not depend on the way in which it 
is described. Let us write the expression of that Gibbs energy in both 
descriptions, and equal those expressions, which gives us: 


G — ^aMa T ^bMb ~ n A A a"* - n A 2 A a, T n B A B [2.42] 

If we look at relation [2.38], we can deduce: 

li^=p\ [2.43] 

and in light of relations [2.37] and [2.39], we deduce: 

M\=Ma [2-44] 

R 'a, = 2/U [2.45] 
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2.5.3. Calculating the extent of the equilibrium in an associated solution 

The law of mass action, written for equilibrium [2.R1] (with a solution 
considered to be perfect) gives us: 

K = = y /2 ,, M-g^) + » B ] P-46] 

(x' A ) n- A (l-a) 

However, the equilibrium constant can be expressed on the basis of the 
partition functions, according to relation [A.3.50] from Appendix 3, which 
gives us: 


K = 


"(A 2 ) 


( Z '»(A)) 


exp 


Ah\T) 
R T 


[2.47] 


In this relation, A r h°(T) is the fundamental energy level of vibration of 
the molecule A 2 . 

Thus, by combining expressions [2.46] and [2.47], we can determine a. 


2.5.4. Calculating the activity coefficients in an associated solution 

Although it is not necessary in order to study our solution, we shall 
calculate the activity coefficients for the components in the first description. 

If we use reference (I) - the pure substance reference - we can write: 

>u A =gl + RT]n/fx A [2.48] 

and 

M\ = g° A +RT\nx\ [2.49] 

Using expression [2.47], we obtain 


[2.50] 
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In the same way, for component B, we find: 

/b )= — [2-51] 

x B 

Using relations [2.35], [2.36], [2.41] and [2.50], equation [2.50] becomes: 

[2.52] 


.//)_ («A+«B X 1 -^) 
I A ~ 


n A (l- a / 2) + n B 
Similarly, for B, we find: 

(«a +«b) 


/b° = 


n A (1 - a / 2) + n B 


[2.53] 


Thus, the two activity coefficients we are seeking can be expressed as a 
function of the extent of the equilibrium [2.R1]. 


2 . 5 . 5 . Definition of a regular solution 

Consider a solution with N components. It is said to be regular if its 
molar mixing entropy is identical to that of a perfect solution, given by 
expression [A.2.10] (see Table A.2.1 in Appendix 2). 

By comparing this expression with the general formula for the molar 
entropy of mixing of a solution: 


S mix = R& 


’ -biy^-lnx. 
3 T 1 


[2.54] 


In view of expression [A.2.10], regardless of the composition, we obtain: 

[2.55] 


+ ln^) = 0 


3 T 


In order for this equation to be satisfied, we must have the following for 
any component v. 


x i(T dlny * + In '/ i I) = 0 
' 3 T ' 


[2.56] 
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Thus, by integrating: 


Tln^' 71 =%(x) [2.57] 

The value thus defined, 9t, (x) , is thus temperature-independent (we have 
discounted the influence of the pressure on liquid solutions). It depends only 
on the composition x. 

Based on the definition we have just given, it is obvious that the excess 
entropy of a regular solution is null, as are the excess partial molar entropies 
of each it its components - i.e.: 

S xs = 0 and S™ = 0 [2.58] 

As the entropy of mixing of a regular solution is identical to that of a 
perfect solution, we can deduce that the molar heat capacity obeys Kopp’s 
law of additivity (see expression A.2.14 in Table A.2.2, Appendix 2). 

We can see that the regular solution model does not give us the variations 
of the activity coefficient with the amounts of substance. In fact, it is a 
family of models which includes several possible subsets - each one defined 
by a relation 9^ (x) . We shall come across two of these subsets in our study: 
the Van Laar equation (see Table 3.3) and the Hildebrand- strictly-regular 
solution model, which we shall touch on. 


2 . 5 . 6 . Strictly-regular solutions 

A group of solutions known as strictly-regular solutions is of significant 
interest for a number of reasons: 

- it yields the most mathematically simple form of a non-perfect and non- 
dilute ideal solution; 

- it demonstrates that complete knowledge of the solution includes all the 
phenomena that are typically supeiposed on the model artificially. Hence, in 
certain cases, it includes the phenomenon of demixing, and the existence of 
azeotropic mixtures; 

- in practical terms, it is applicable in numerous cases in metallurgy with 
mixtures of liquid metals, and in chemistry with solutions of homologous 
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compounds. It is found in the modeling of other phenomena such as the 
concept of order/disorder in alloys, or the isotherm of non-“Langmuirian” 
adsorption of a gas on a solid; 

- we have a microscopic model of strictly-regular solutions, as we do for 
perfect solutions, which enables us to demonstrate an initial cause of 
imperfection of the solutions. 


2.5.7. Macroscopic modeling of strictly-regular binary solutions 

Let us take the Margules expansion, limited to the first non-zero term. In 

light of relations [2.2] to [2.5], assuming that /?, = /? 2 = — , where B is a 

parameter independent of the temperature and of the composition, the 
activity coefficients will be written as: 

foyl I)= j x 2 and M'’ = j x f [2-59] 

In view of the definition of B, it is clear, given [2.57], that such a solution 
is regular. 

Note that the function of the pressure is not specified, and therefore the 
model includes all solutions defined by these two equations, whatever 
the function B(P) and, in particular, B can very well be independent of the 
pressure. 

The solutions that we have just defined constitute a family of regular 
solutions, known as strictly-regular solutions in the sense of Hildebrand. 

In these solutions, the two components play a perfectly symmetrical role. 

As strictly-regular solutions are, by definition, regular, we can deduce 
that the entropy of mixing is the same as that of a perfect solution and 
therefore that the excess entropy is null. 

We shall now calculate the other main thermodynamics involved in 
strictly-regular solutions. 
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The excess partial molar Gibbs energy of a component is obtained 
directly on the basis of the activity coefficient: 

G xs = R T In y = R Bx\ and G xs = R T In y, = R Bxf [2.60] 

The total molar excess Gibbs energy will therefore be: 

G xs =x,G^ + x 2 GT = RRxjX 2 [2.61] 

As the excess entropy is null by definition, the excess enthalpy is equal to 
the excess Gibbs energy, and thus for the partial molar values, we would 
have: 


H; s = G xs = R Bx\ and /7 XS = G xs = RBxf [2.62] 

and for the molar excess enthalpy: 

H xs =xjif + xj{f = R Bxyx 2 [2.63] 

As the enthalpy of mixing for a perfect solution is zero, the enthalpy of 
mixing for a strictly-regular solution is identical to its excess enthalpy, so for 
the the partial molar values: 

= RBxl and H?* =H™ =RBxf [2.64] 

and for the molar enthalpy of mixing: 

//;; x = x, ~h^ + x 2 Trr = r b Xi x 2 [2.65] 

We can see that the plot of the mixing enthalpy is a parabola (see 
Figure 2.2) as a function of the molar fraction of a component because: 

TT mix 

H m 


R B 


= Xi(l-Xi) 


[ 2 . 66 ] 
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r mix 
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Figure 2.2. Parabolic form of mixing enthalpy 
for a strictly-regular solution 


The Gibbs energy can immediately be deduced from the definition: 



[2.67] 


In the knowledge that the mixing enthalpy is the same as that for a perfect 
solution, we find the following for the molar Gibbs energy of mixing: 


G,™ x = R Bx x x 2 + RT(x t In x, + x, lnx 2 ) 


[ 2 . 68 ] 


and for the partial molar values: 


G, mix = R Bx] + R Txj In x, and G 2 m,x = RBxf + R Tx 2 In x 2 


[2.69] 


As we saw earlier, because a strictly-regular solution is a regular solution, 
it obey’s Kopp’s law of additivity (see relation A.2.14 in Table A.2.1, 
Appendix 1) for the molar heat capacities. 

We can easily calculate the expression of the activity of a component at 
temperature T. Immediately, we find: 



\ 


fl 2° = = X 2 eX P 


y 


[2.70] 
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The relation is absolutely symmetrical for component 1. 



Figure 2.2. Activity (in convention I) of a component 
of a strictly regular solution as as a function of the composition 


Figure 2.3 illustrates the activity of a strictly-regular solution as a 
function of the composition. We can see that the solution displays positive 
deviation in relation to Raoult’s law if B/T < 0, and negative deviation if 
B/T > 0. If B/T is greater than 2, the inflection point indicates demixing into 
two solutions, whose compositions are given by the endpoints of the plateau. 
The activities on either side of the plateau correspond to a metastable phase. 

The curves for the second component are, obviously, symmetrical. 

NOTE 2.2.- the strictly-regular solution becomes a perfect solution if the 
parameter B is 0 for all compositions (B/T = 0). The corresponding curve in 
Figure 2.3 is Raoult’s straight line. 


2.5.8. Extension of the model of a strictly-regular solution to solutions with 
more than two components 

Certain authors prefer to choose relation [2.89] from the expression of the 
excess Gibbs energy as the definition for strictly-regular solutions. This 
second definition is rigorously identical to that which we have chosen 
(relation [2.87]), because relation [A.2.32] (see Table A.2.2 in Appendix 2) 
li nks the activity coefficients and the partial molar Gibbs energy values. 
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On the basis of this definition, we can extend the concept of a strictly- 
regular solution to systems with more than two components, but preserve the 
symmetry that is present in binary systems. For instance, for a system with N 
components, we can say that a solution is strictly regular if its excess molar 
Gibbs energy takes the form: 

G ”=iiz p- 71 ] 

;=1 j = i + 1 

The coefficients Ay are positive coefficients, independent of the 
temperature and composition. 

By application of this relation to a two-component system, obviously, we 
obtain relation [2.60], with: 

4 = RB [2.72] 

For a system with three components 1, 2 and 3, we find: 

t 'j A^^XyX^ I I [2 . V3 ] 

Based on the derivative of this function in relation to one of the quantities 
of matter, we immediately obtain the activity coefficient of that component. 
For example, for the three-component solution, we obtain: 

RTln y[ I) = = RT(4,jc, + 43 * 3 ) [2-74] 

OTly 

Multicomponent strictly-regular solutions constitute an excellent model 
for mixtures of molecules of the same polymer in different states of 
polymerization, which is understandable in view of the similitude of the 
molecules involved. 


2.6. Athermic solutions 

Alongside regular solutions (see section 2.5.1), we can define athermic 
solutions, which are primarily found when amorphous or molten polymers 
are dissolved in solvents. 
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2 . 6 . 1 . Thermodynamic definition of an athermic solution 

In the same way as we defined a regular solution as one which has the 
same entropy of mixing as a perfect solution, we shall define an athermic 
solution as one which has the same enthalpy of mixing as a perfect 
solution - i.e. zero. Of course, its excess molar enthalpy is also null. 


2 . 6 . 2 . Variation of the activity coefficients with temperature in an athermic 
solution 


By applying the above property of definition to expressions [A.2.24] and 
[A.2.8] for a perfect solution, we deduce: 


X 


d In 
dT 


= 0 


[2.75] 


As the activity coefficients do not depend on the temperature, they will be 
functions only of the composition and perhaps of the pressure. Therefore, 
this type of solutions is referred to as athermic. 


ff = f‘\p, x ) [2.76] 

Instantly, we can derive the following equality from this: 

M=l^ [2.77] 

Thus, the partial molar enthalpy of each component is identical to its 
molar enthalpy when it is pure (in the same state of aggregation). 


2 . 6 . 3 . Molar entropy and Gibbs energy of mixing for an athermic solution 

Using the general expression of the excess entropy [A.2.33], for the 
excess molar entropy of an athermic solution, we obtain: 


[2.78] 
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Hence, the molar entropy of mixing is: 


Sr=-R$>,( In^ + ln*,-) 


[2.79] 


As the molar enthalpy of mixing is null, the molar Gibbs energy of 
mixing is given by: 


G mix 
m 


- ts : 


R/^ln^Tln*,) 

i=l 


[2.80] 


2.6.4. Molar heat capacity of an athermic solution 

The excess molar heat capacity is written simply as: 

r)H xs 

C X P S = — — = 0 [2.81] 

dT 

Thus, the molar heat capacity is the same as that for a perfect solution and 
therefore obeys Kopp’s additivity law [A.2.14]. This law can be said to be 
fairly widely applicable, because we have seen it at work with perfect 
solutions, regular solutions, athermic solutions and approximately with 
dilute ideal solutions. 

NOTE 2.3.- a solution which is both regular and athermic is a perfect 
solution, because these two solutions have the same enthalpy and the same 
entropy, and therefore the same Gibbs energy. Hence, they are identical. 

In this chapter, we have examined a number of simple analytical models 
which are extremely useful for the study of equilibrium states. However, 
these models involve parameters with unknown physical meaning. In 
addition, they are often not sufficiently accurate. For that reason, it is worth 
considering the modeling of solutions by a microscopic description. 
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Microscopic Modeling of Liquid 
Molecular Solutions 


The purpose of modeling liquid solutions is to understand the structure 
and the properties of the medium, and to help predict certain behaviors 
by giving expressions for the molar Gibbs energies or activity 
coefficients. These data can then be used in practical applications such as 
the establishing of chemical equilibrium states or the plotting of phase 
diagrams. 

Solutions are complex media, and usually the model needs to be 
“calibrated” on experiments in order to determine certain constants. The 
advantage to using a model, in terms of the predictive aspect, is obviously 
that it requires a minimal number of terms to be determined, and therefore a 
minimal number of experiments, to calibrate it. Nevertheless, this number 
must be sufficiently high for the model to represent reality as closely as 
possible. 

In the existing body of literature, a great many different microscopic 
models are put forward. Hence, it is impossible to discuss all of them here. 
Instead, we shall limit ourselves to a presentation of the most widely-used 
models, which are practical and illustrate modem ideas by way of the 
approaches and hypotheses that they implement. Certain models will not be 
examined in detail because, whilst they were groundbreaking when they 
were first published, today they are known to be specific cases of other, 
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more general, models. We shall content ourselves with listing these models 
in Table 3.3. 


3.1. Models of binary solutions with molecules of similar dimensions 

Suppose we have a binary liquid solution - a mixture of two components 
A and B, which contains n A moles of component A and n B moles of 
component B, which equates to N A molecules of A and N B molecules of B. 
This solution satisfies the following six conditions: 

1) There is an order over a short distance, and the short-distance lattice is 
the same for the solution as for the two liquids A and B. Thus, these three 
liquids exhibit the same coordination number z (e.g. 12 for the hexagonal or 
cubic lattice with centered faces). 

2) The molar volumes v° and v° of the two pure liquids are sufficiently 
close. Certain authors have evaluated that this condition was acceptable if 
the ratio of the radii of the two molecules A and B, which are supposed to be 
spherical, was no greater than 1.25; 

3) The free volumes v /(A) and v f(B) of the two pure liquids (see section 
1.3.1) are near to one another (to within 30%); 

4) The molar volumes and the free volumes of the two liquids are 
not altered by the operation of mixing, and therefore the mixing volume is 
null, which means that the volume of the solution is given by the additive 
law: 


^ = «A V A+«B V B [ 3 - 1 ] 

5) The potential energy of interaction Ej can be considered to be the sum 
of the contributions of the closest neighboring pairs of molecules; 

6) The model of the liquid chosen is Guggenheim’s smoothed potential 
model (see section 1.3.1), both for the pure liquids and for the solution, and 
we use £ A{A) and e B(B) to represent the values of the lattice energy £ 

which appears in expression [1.22], respectively for each of the pure liquids 
A and B. 
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The mean energy of interaction of two molecules of A is -2f A(A) / z , and 

similarly for two molecules of B it is ~2e B(B) / z . Hence, we can say that the 

energy levels of pairs of molecules A-A and pairs de molecules B-B are, 
respectively: 

^-^and^- 2 ^ [3.2] 

z z 

NOTE 3.1- As we are essentially interested in the mixing values, the pure 
liquid model chosen has no impact on the results. 

Let us define an energy value, called the exchange energy > w AB , such that, 
if we begin with the two pure liquids A and B, the exchange of a molecule of 
A and a molecule of B between the two liquids increases the energy of the 
system by the amount 2iv ab . During this process, z A-A pairs and z B-B 
pairs have been destroyed, and 2 z A-B pairs have been created. 

Thus, the mean potential energy of the pairs A-B will be: 

£ ab ~ (~2£ A ( A ) — 2£ B ( B j + 2w ab ) / 2z = £ A ( A j - £ B ( B ) + w AB j/ z [3.3] 


In view of equations [3.2] and [3.3], the exchange energy is: 


w 4n = z 


( p i p \ 

^ C AA ' C BB 

^4R 


[3.4] 


The property vv A u will be temperature-independent, so long as the type of 
environment of a molecule and therefore z are not altered. 

NOTE 3.2.- The energy of interaction between a molecule of A and a 
molecule of B in the solution can be written as £ m = — 2£ AB(AH) / z. Thus, 

the exchange energy can be written as follows, if we take account of 
relation [3.3]: 

W AB — ^A(A) ^B(B) _ ^AB(AB) 



64 Modeling of Liquid Phases 


We can establish a general form of the configuration integral in the case 
of interactions by pairs. 

If we look again at relations [A.3.25] and [A.3.37], we can write the 
following formula for the translational partition function of a substance A: 
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[3.5] 


where co A is the set of position coordinates of the N A molecules of A in a 
given configuration. 

The integral of equation [3.5] is thus extended to all configurations of the 
N a molecules de A. 


Similarly, for a binary mixture of two species A and B, we can write the 
canonical partition function of translation in the form: 
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[3.6] 


This time, the integral is extended to all configurations of Na molecules 
of A and N B molecules of B (the extension to solutions with more than two 
components is obvious). 

We shall use Z A ^ ;/ . to denote the canonical partition function of the pure 

component A without interaction, i.e. the perfect gas; similarly, 

denotes that of component B. The canonical partition function of translation 
of the solution can be written as: 
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[3.7] 
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The problem of modeling of the solution involves being able to express 
the configuration integral of the mixture, whose value is: 


4b = 


l l 
7VJAV 


Ml 


exp 


,j4_ 

v k B T J 


d <y" A d<yg B 


[3.8] 


For each configuration, the solution comprises a certain number of pairs 
A-A, pairs B-B and mixed pairs A-B. 


Let us look at a particular configuration of the solution in which the 
number of mixed pairs A-B is zX. 


Thus, the number of neighbors of a molecule of A which are not 
molecules of B is z( N A - X) . Hence, the number of pairs A-A 
is z(jV a - X) / 2 and, by the same token, the number of pairs 
B-B is z(N B - X) / 2. Therefore, the total number of pairs N p of the three 
species is: 


v, =| (w„ +N„ + V„) = £(V a + AT.) 


[3.9] 


For this particular configuration, characterized by the value of X, the 
potential energy of the solution, due to the interactions, will be: 


4 = 


z(N a -X) 


2e 


A(A) 
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z(N b -X) 


2f, 


B(B) 


+ zX- 


. ( 4(A) 4(B) + 4b) 


[3.10] 


Thus: 


4 — 4.4(a) 4 £ b(b) + Xw AB 


[3.11] 


This relation can also be written as follows, taking account of 
expressions [3.2] and [3.8]: 


4 = 


zN ,£. 


- + - 


zN a £ D 


-+Xz 


[3.12] 
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Note in passing that, if x A and x B are, respectively, the molar fractions of 
the components A and B, the internal energy per mole of solution 
(iVA+iVB=N a ), which is also its enthalpy, because the volume of mixing is 
null, is: 


U m =H m = N a 


ZX A&AA | ZX B^BB 


+ Xz 


^AA ^BB 


[3.13] 


By substituting expression [3.12] back into equation [3.8], the exchange 
integral becomes: 
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In order to calculate the integral of expression [3.14], we define the 
value Y by relation: 
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As, according to hypothesis 3, the free volumes v f(A) and v /ui) are very 

close, we can replace the arithmetic mean with the geometric mean, and use 
the following approximation: 


N A V f(A) + N 2 V AB) _ / N. N B \{ N a +N b 
= VnA) v f( B) 11 


n a + n b 




[3.16] 
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In light of relation [3.16], we can integrate the left-hand side of 
expression [3.15], and we obtain: 
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exp 


Yw 


V k s T J 


The configuration integral / AB can be broken down into three terms: 
an interaction integral for the pure liquid A (term given by relation [1.22]), 
an interaction integral for the pure liquid B and a mixing integral 
which is: 


t mix 
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[3.18] 


The composition Y appears as the fictitious number of mixed pairs that 
would yield the right value of the partition function if the molecules were 
distributed at random. 


Note 3.3 - As pointed out earlier, the choice of the model of a pure liquid is 
no longer relevant in the expression of the integral of mixing. 

The integral of mixing is the contribution of the mixing to the system’s 
overall partition function. 

We can now calculate the Helmholtz energy of mixing and the Gibbs 
energy of mixing, which are equal because the volume of mixing is null, 
with the following formula: 
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[3.19] 
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Using Stirling’s approximation, this Gibbs energy of mixing takes the 
form: 


(jmix _ pmix _ y w 


-k B r 


N a In 
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[3.20] 


Thus, we need to calculate the value Y. In order to do so, we write (X) 
for the value of X at equilibrium, which is such that: 
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[3.21] 


By deriving equation [3.15] in relation to temperature, we can verify the 
expression: 


(. X) = Y-T 


dY 

dT 


d [Y/T] 
0(1 IT) 


[3.22] 


Thus, ultimately, calculation of the partition function of our model of a 
solution boils down to determining (X) and then Y using expression [3.22]. 

Finally, by feeding the value Y thus obtained back into expression [3.20], we 
find the value of the Gibbs energy of mixing. 


3 . 1 . 1 . The microscopic model of a perfect solution 

We have already given (see relation [A.2.28] and Table A.2.1 in 
Appendix 2) the macroscopic representation of a perfect solution and the 
resulting properties. Now, we are going to construct a microscopic model, to 
see which hypotheses yield the same expression for the Gibbs energy of 
mixing given by expression [A.2.6] (see Table A.2.1). 
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We return to the previous model (section 3.1), and we shall now 
determine Y by adding a seventh hypothesis to our model: the exchange 
energy is null, i.e.: 


2f — F 4 - F 

... _ Z ' C AB e AA ' C BB _ A 

W A B U 


[3.23] 


Thus, by substituting this value into relation [3.13], we see that the total 
enthalpy is equal to the sum of the enthalpies of the two pure liquids, and 
therefore the enthalpy of mixing is zero. This means that our two solutions 
are mixed with no alteration of the energy - i.e. with no thermal effect. 


With our new hypothesis, relation [3.18] becomes: 


rmix 
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[3.24] 


The Helmholtz energy of mixing is the same as the Gibbs 
energy of mixing (as the mixing volume is null), and in light of 
relation [3.19] (w^ = 0) , we calculate: 


G mix = k B T 
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[3.25] 


In order to find the molar value, we make /V A + AT = Na and n A + n B = 1, 
so that x A = « A an d x B = n B . Thus, the molar Gibbs energy of mixing is: 

G'T — F,T = RT ( x A In x A + x B In x B ) [3.26] 


We can see that this relation is indeed identical to expression [A .6] (see 
Table A.2.1 in Appendix 2). Thus, we have here a model of the perfect 
solution, and note that this solution is not, as might be imagined, a solution 
in which there are no inter-molecular interactions, as is the case with a 
perfect gas. Rather, it is a solution in which the energy of interaction f AB 
between two molecules of A and of B is the arithmetic mean of the energies 
of the A-A and B-B pairs because, in view of relations [3.4] and [3.23], we 
obtain: 


p _ ^AA ^~BB 
fc AB 


2 


[3.27] 
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Of course, as is often the case in modeling, we have established a model 
of a perfect solution, but there is nothing to say that this is the only possible 
model. 


3.1.2. Microscopic description of strictly-regular solutions 

We shall now return to our model of a solution given in section 3.1, with 
the six hypotheses upon which it is founded, but exclude hypothesis 7, given 
by equation [3.23]. The exchange energy w A n is no longer null. 

NOTE 3.4.- Whilst we continue to state the hypothesis of a random 
distribution of the molecules of A and B, this hypothesis appears to 
contradict the existence of an exchange energy, whose presence means that 
the minimum of the Helmholtz energy cannot correspond exactly to the 
random distribution. We get around the contradiction by accepting that the 
short-distance order that would be caused by this exchange energy is 
annihilated by the thermal agitation. In other words, we accept the condition: 

W AB« k B 7 ’ [3.28] 

Thus, let us write the random distribution of the mixed pairs A-B. Around 
a molecule of A, there are, on average, zx B molecules of B, so the mean 
number of such pairs will be given by: 

z { x ) = N A ZX H = N b zx a [3.29] 

This enables us to write: 


(X ) 2 = N a N b x a x b =N a N b 
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[3.30] 


Thus: 

WM^aa-WX^bb-W) [3.31] 

Thus, is independent of the temperature. This is what is known as 
the zero-order approximation or the Bragg-Williams approximation. 
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This temperature-independence of means that, in light of relation 
[3.22], we can write: 

(X) = Y [3.32] 


and, using relation [3.29], we deduce from this that: 


(*) = * 


N a N b 

n a + n b 


[3.33] 


The integral of mixing is obtained from relation [3.18]. It is of the same 
form as relation [3.18]. It will be the product of two terms: one identical to 
that in relation [3.24], pertaining to the perfect solution, and the other 
comprising the exponential, which will, in fact, be the excess term. 
According to equation [3.33], this term is: 


G" = Yw. 


N a»A»B 
»A + »B 


W 4 


[3.34] 


and for the corresponding molar value («a+«b = 1), we obtain: 


Gm = N a*A X B W AB 


[3.35] 


As Y is independent of the temperature, the excess entropy is null, and 
we find the following molar enthalpy: 


H” =~T 


3(G-/r) 


dT 


= N ,Vb»'ab 


[3.36] 


Thus, we see all the properties of the strictly-regular solution 
(section 2.5.6) and by identification between equations [3.35] and [2.61], we 
can deduce: 


B = 


N w A 


R 


[3.37] 


The difference between a strictly-regular solution and a perfect solution 
lies in the exchange energy, which is null in the case of the latter. 
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Note 3.5 - the null excess entropy is due to the temperature-independence 
of Y, which means, in view of relation [3.22], that the equilibrium value 

is also independent of T, and we have shown that the random 

distribution of the molecules satisfies this condition. There is nothing to 
show that only this distribution leads to a null excess entropy, which 
explains why not all regular solutions are strictly regular. 

Certain authors work differently to construct the model. They accept the 
hypothesis of random distribution of the molecules. They then deduce that 
the excess entropy is null. Knowing the excess enthalpy, given by relation 
[3.13], they then calculate the Gibbs energy, applying relation [3.36] in 
reverse. Integration is then performed, making the hypothesis that, if the 
temperature tends toward infinity, the mixture behaves like a perfect 
solution, which is to say that its excess enthalpy is null. In fact, this way of 
working seems much simpler, and it does yield the right result. However, the 
assertion that the mixture behaves like a perfect solution if the temperature 
tends toward infinity is not easily conceived of. Thus, we prefer the previous 
construction. This presentation also has the advantage of being able to be 
used in the case of another model (see section 3.2.4). 

3.1.3. Microscopic modeling of an ideal dilute solution 

Let us look again at the above model of a strictly-regular solution 
(section 3.1.2). If the number of molecules of B is small in relation to that of 
A - in other words, if the solution is very dilute in terms of B content - then 
relation [3.33] gives us: 

(X) = Y = N B [3.38] 

Relation [3.34] yields: 

G “= N a«B W AB [3-39] 

The activity coefficient of B in convention (I) thus satisfies the equality: 

'Ty'-rXS 
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[3.40] 
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Therefore, we can deduce that: 


Zb’ = exp- 


N a^AB 

RT 


[3.41] 


This activity coefficient is independent of the composition of the 
solution; it depends only on the temperature. 

The excess molar Gibbs energy is thus: 

G“ = N a w AB =-Rr« B lny'" [3.42] 

In this excess Gibbs energy, the solvent makes no contribution at all. This 
tells us that y A ’ = 1 . Consequently, as with the solvent, conventions (1) and 
(II) are identical. We deduce that: 

Ti" = rT = 1 [3-43] 


In convention (II), therefore, the activity coefficient of the solute B is also 
constant, and as with the infmitely-dilute solution, the value of this constant 
is 1. From this, we deduce that: 

= 1 [3-44] 

Here, we encounter all the characteristics of the ideal dilute solution 
studied in section 2.4. 


Henry’s constant is therefore given, on the basis of relation [2.23], by: 


k bh = 7b ] = exp 


N a^A] 

R T 


[3.45] 


Thus, the ideal dilute solution is the particular case of the strictly-regular 
solution wherein there is a very small proportion of one of the components. 

From expression [A.2.21], we deduce the physical meaning of the 
chemical potential of the solute for the reference state (II): 


M7 = S°s - N a w AB 


[3.46] 
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The chemical potential of reference state (II) is that of reference state (I) 
less the molar exchange energy. 

The most general of the three models of solutions that we have just 
examined is that of the strictly-regular solution. The other two models seem 
to be two particular cases; that of the perfect solution if we make w AB = 0 
and that of the ideal dilute solution if we make n A » n B . In all these 
models, we assumed a random distribution of the molecules, which is a 
fairly restrictive hypothesis, as we mentioned in our discussion of the 
strictly-regular solution. It is this hypothesis which we are going to examine 
in greater detail in the next section. 


3.2. The concept of local composition 

In the model of strictly-regular solutions, we have indicated a 
contradiction in the fact that the A-B bond involved a different amount 
of energy than the A-A and B-B bonds, but admitted that this difference 
did not prevent random distribution of the pairs A-B - in other words, 
the local composition of the liquid around a molecule of A or a molecule 
of B was the same as the overall composition. Now, we shall refine 
our models by examining local distributions different to the overall 
composition. 


3.2.1. The concept of local composition in a solution 

In reality, the exchange energy has consequences for the local distribution 
of the molecules. Indeed, it is easy to see that if there is a stronger force of 
attraction between a molecule of A and a molecule of B than there is 
between two molecules of A or between two molecules of B, then a 
molecule of A will tend to be surrounded by molecules of B, and therefore 
the composition of the immediate environment of the molecule of A will not 
be the same as the overall composition of the solution - it will be richer in 
molecules of B. The opposite effect would be observed if the two molecules 
A and B had a lesser force of attraction than two molecules of A and two 
molecules of B. This is the concept of local composition, first expressed by 
Wilson. 
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In order to describe this concept of local composition, we use a 
nomenclature that enables us to clearly distinguish, every time, between the 
central molecule and its closest neighboring molecules. 

Ny denotes the number of molecules of i around a molecule of j, and we 
define the conditional probability Py that a molecule of i will be located in 
the environs of a molecule of j by the expression (where N denotes the 
number of components in the solution): 

N, 

P v =^~ [3-47] 

IX 

1=1 

Note that this probability also represents the local molar fraction. 

Evidently, we have: 

IX =1 [ 3 - 48 ] 

1=1 

Thus, for a binary solution, we distinguish: 

- P AA : local molar fraction of the molecules of A around a molecule 
of A; 

- P BA : local molar fraction of the molecules of B around a molecule 
of A; 

- P BB : local molar fraction of the molecules of B around a molecule 
of B; 

- P AB : local molar fraction des molecules of A around a molecule of B; 
with the equalities: 

P aa + P ba = 1 [3.49a] 

and P BB +P AB = 1 [3.49b] 

We shall suppose that it is always possible to define a weighting factor ky 
such that the local composition is the overall composition weighted by that 
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factor. In other words, we shall define the factors k tJ by way of the following 
relations - e.g. for a binary solution (x A and x B are the overall molar fractions 


two components in the solution): 


/k 

~—k 

a BA 

[3.50a] 

^AA 

X A 


P\n . 

= —k 

ft 'AB 

[3.50b] 

P BB 

x B 



By combining equations [3.49a], [3.49b], [3.50a] and [3.50b], we obtain: 

[3.51] 


P AA = 


X A X B^BA 


and by substituting the result back into expression [3.50a], we find: 

x D k a . 


Pba = 


X A X B^BA 


[3.52] 


Similarly, around molecule B, we find: 


/> =■ 


X B X A^AB 


and: 


_ X A^AB 

MR ' 


X B X A^AB 


[3.53] 


[3.54] 


NOTE 3.6.- If the arrangement of the molecules of A and B is random, 
obviously we have P AB = x A and P BA = x B , and therefore k m = k BA = 1 . 


3.2.2. Energy balance of the mixture 

We shall now recalculate the energy of mixing of /V A molecules of A and 
N b molecules of B in that context of local composition, by finding the 
energy balance of that mixture. We suppose that in pure liquid A, a molecule 
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of A has a coordination number z A , meaning that it is surrounded by z A 
molecules, which exert an influence on it in energetic terms. Similarly, in 
pure liquid B, the coordination number is z B . In solution, we use z \ and z’ B 
respectively to denote the coordination numbers of the molecules of A and 
the molecules of B in the solution. There is no reason, a priori, for all these 
coordination numbers to be equal, which runs counter to the hypothesis we 
have employed up until now. 

The interaction between two molecules that are not immediate neighbors 
is ignored. 

In order to extract a molecule of A from the pure liquid A, we must break 
z A A- A bonds and therefore provide the energy z A Aa- We place that 
molecule in the pure liquid B. This transfer forms new bonds and therefore 
releases the energy z A (v A aAa + TbaAa)- The variation in energy due to the 
transfer of a molecule of A from the pure liquid A into liquid B will 
therefore be: z A e AA - z ’ A (y AA e AA + y BA £ba)- Obviously, we can calculate a 
symmetrical expression for the transfer of a molecule of B from liquid B into 
liquid A. The variation in internal energy to form a mole of mixture, 
therefore, will be: 


jj-mix 


Na [a [ Z A A\ z a(TaaAa “*“TbaAa)J 
] 

^ +A [ z bAb — z b (TbbAb + TabAb )] 


[3.55] 


Taking account of relations [3.50a], [3.50b], [3.53] and [3.54], we find 
the following for the excess- or mixing molar internal energy: 


, jxs _ Na AA^ba z (g BA A x b^ab z b(Ab Ab) [3 561 

m s-t T7- ts- L • J 


A+A^ba 


A+A^ab 


Based on our knowledge of the internal mixing energy, we calculate the 
Helmholtz energy using the relation: 


9(i/r) 


[3.57] 
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Integration gives us: 


FI 


F x 


u: d 




[3.58] 


If T lim is infinite, at this temperature F/T m becomes independent of the 
temperature because the mixture is identical to a mixture of perfect gases 
without interaction between molecules, but FIT m may depend on the 
composition or the density of the liquid. 

Finally, the modeling of the solution, i.e. the establishing of the excess 
molar Helmholtz energy function, involves expressing the value of the 
coefficients ky from relation [3.56], and most of the existing models differ 
only in this evaluation. 


3.2.3. Warren and Cowley’s order parameter 

For a solution containing molecules i and j, Warren and Cowley defined 
an order parameter TJ^ as being: 


nb of molecules of i around a molecule of j 
mean nb of molecules of i around a molecule of j 


In the case of a binary mixture, therefore, we would have: 





with the condition: 


0</7ab=£1 

and in parallel: 

77 ba = 1-— and 0 </ 7 ba <1 

x B 


[3.59] 


[3.60] 


[3.61] 


We can obviously deduce that: 



^BA — *b(1 %a) 

an d ^AB — X A ( 1 _ VaB ) 
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[3.62a] 

[3.62b] 


NOTE 3.7.- if the distribution of the molecules is random, then P AB = x A and 

%A=1- 


The two order parameters // Al! and // liA are connected; indeed, the 
number of pairs A-B can be expressed in two ways, such that: 


N N =N N 

JV AB JV B 


Hence: 


z' P N =z' P N 

* A 1 BA JV A B 1 AB JV B 


[3.63] 


[3.64] 


By introducing the values of P BA and P AB in equations [3.62a] and 
[3.62b], we find: 


Z A^a(! 7ba) Z B X aNr(\ Vab) 
If we note that we have the equality: 

n a n» 


X &N a — — 


n a + n b 


[3.65] 


[3.66] 


We obtain: 


Z A %A Z B >7aB Z B Z B 


[3.67] 


Hence, the two parameters of opposite orders are not independent. 

If the coordination numbers of two types of molecules are the same in the 
solution ( z' A = z' B ) , then the two Warren-Cowley order parameters are also 
equal ( r/ AB = tj BA ) . 


There is an obvious relation between the order parameters and the 
weighting factors k,f, indeed, according to equations [3.51], [3.52], [3.62a] 
and [3.62b], we can write: 
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n ba = 1 ^ BA . an d ?1 AB = 1 [3.68b] 

•*A X B^BA X A X AB X B 

This expression shows that if =0 or & AH = 1 , the distribution of 
molecules A and B is random. 

If // Al! > 0 or k m < 1 , this means that P AH < x A . Thus, on average, there 
are fewer mixed pairs than in the disordered solution, so there is a tendency 
toward dissociation of the pairs A-B. On the contrary, if rf m < 0 or k m > 1 , 
that would mean that there is a certain tendency toward the association of 
molecules of A with molecules of B. 

3.2.4. Model of Fowler & Guggenheim ’s quasi-chemical solution 

Let us return to our binary solution formed of two components A and B. 
With the exchange energy being given by vv A b, we can write that the 
system’s partition function of exchange is the sum of all the values of X, 
with each energy state being weighted by a degeneration g(X ) , so that: 

x 

g(X) is the number of arrangements which are possible for a given level 
of energy. For a given value of X, the different pairs of atoms may be 
arranged in accordance with the number of possibilities, proportional to: 

(zN/2)\ 

[z(Af A -X)/ 2 ] ![ z (tf B -*)/2] ![zAab/2] ![zAW2] ! 

Thus, we can write: 


(2zXw m zN^S^ zN m £ BBi 
zk B r 


[3.69] 


g(X) = k 


(zN/2)\ 

[z(N a -X)/2]![ z (A b -X)/2] ![zAab/2] ![zAba/2]! 


[3.70] 


In this expression, the first and second term in the denominator, 
respectively, are relative to the numbers of pairs A-A and B-B, and the third 
and fourth term respectively relate to the numbers of pairs A-B and B-A. 
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We have introduced a coefficient of proportionality, because the A-B 
bonds are not totally independent of one another. Figure 3.1 shows an 
example of this non-independence. Indeed, suppose that on the lattice (which 
we shall take to be flat for the purpose of simplicity), we have placed the 
molecules A a , B b , and B c respectively on sites a, b and c. If we place a 
molecule B on site d, it is clear that the bond between sites a and d is 
necessarily a mixed bond A a B d . In other words, the pairs are not distributed 
entirely randomly. The coefficient k depends on the lattice chosen to 
represent the cell. 



Suppose we use g a ( X ) to denote the possible number of arrangements of 
pairs when they are distributed “at random”. According to equation [3.70], 
we have: 




[z(^a-A*)/ 2 ] 


(zN/ 2) ! 

l[z(A b-A*)/2] ![zAab/ 2 ] ![zAba/2 ] ! 


[3.71] 


However, we know that, for a random distribution, this number is: 


g a (*) = 


N\ 

K] M • 


[3.72] 


By multiplying and dividing relation [3.70] by the values of g(X) given 
respectively by relations [3.72] and [3.71], which are equal, we obtain: 

ri.r)J -("»- i ")' 2 lh( iy -- r >2l ! r ! n.73, 

[z(^a-T)/2] ![z(^b-T)/2] ![^T/ 2 | i[zX/2] I [N a ] ![IV b ] ! 


To simplify the calculations, we apply Boltzmann’s law (see 
Appendix 3), whereby the sum of all the states in equation [3.69] can be 
replaced by its maximum value. 
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In this case, the number of mixed pairs will be such that: 


^{/ab} 


[3.74] 


(Xj is the value of X which satisfies relation [3.74]. 


We obtain: 


(X) 2 =(A A -(X))(A B -(X))exp-- 


[3.75] 


A laborious calculation performed by Christian in 1975 gives us: 


/ Y~\ — 

V 1 {N A + N B )fi q + 1 


[3.76] 


If we introduce the molar fractions defined by: 

n. N. , n n AC 

x A =- “ — r = 7 — -andx B =- 5 — 7 = 7 — r [3.77] 

(«a+«b) {N a +N b ) («a+«b) (^a+^ b ) 

Thus, for a mole of mixture ( N A + N B - N a ) , we have: 


/ y \ — ^^aTA~B 

' ' N a /?,+l 


[3.78] 


The term [3 q is defined by: 


r \ 

1 + 4 N a N b exp — H ' ar - 1 

L z K T 


= t ( 1 - 2 x a ) 2 + 4x a x b exp- : 


1/2 


[3.79] 
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By substituting expression [3.78] into equation [3.22], we obtain: 


Y = 


n a n b zk B r 2WAB | 

r 2 a\ 

f 2 w ab "] 

N A +N B 2w 12 j 

> 

1 

c 

v zk B T y 


[3.80] 


The lower bound of the integral of expression [3.80] is determined as 
follows: if T — we must observe a random distribution (because then 
w 12 « k B T), in which Y is more or less equal to (2f), according to relation 
[3.22], and where the difference Y-^X'j is finite. Thus, Y IT tends toward 
zero if 1 IT tends toward zero. In view of the definition of /3 q (relation 
[3.79]), we have: 


exp 


"20 

V zk s Tj 




_ rq 


4x a x 2 


Thus: 

d 


^2w AB ^ 


V Zk B T J 


2 A d P q 


(^-1 + 2 x a )(^+ 1 - 2 x a ) 


[3.81] 


[3.82] 


By introducing this into expression [3.80], we find: 


zkj 


n a n b 


2vv. 


3-1 + 2x a 


Xa i n ' q A + Xr In— 2A. _ i n T ? 


A+1-2x a , 3+1 


2x, 


2x a 


[3.83] 


In the knowledge that: 

(N A +N B =\) [3.84] 

By feeding this back into relation [3.20], we obtain the excess molar 
Gibbs energy: 


G xs = —zRT 
2 


P ' +x B -x A 

x n In 


1 (A / + 1 ) 


+ x 4 


, Pq +X A~ X B 

In— L - - 

x A (P q + 1 ) 


[3.85] 
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Note that, as was the case with the strictly-regular solutions model, our 
new model preserves a result that is symmetrical in relation to the two 
components A and B. 

Figure 3.2 shows the difference between the excess molar Gibbs energy 
of the strictly-regular solution and that of this model for the values 
T = 800 K, z = 12 and Naw A B = 30 kJ. The two curves reach their minimum 
point at x A = x B = 0.5. 


G"(Kj/nwle 



Figure 3.2. Comparison of the excess Gibbs energies of a strictly-regular solution and the 
quasi-chemical model (reproduced from [DES 10], p.62 - see Bibliography) 


We can see that the solution is no longer regular. 


The activity coefficients of the two components become the following, by 
derivation of the excess Gibbs energy in relation to «„ remembering that j3 q , 
x A and x B are functions of N A and /V K (the calculation is complicated but 
yields a direct result): 


y(D _ 

/ A 


Pq + *A ~ *B 

*a(A, + 1 ) 


[3.86] 




Pq +X E~ X A 
*B (Pq + !) 


V 


y 


[3.87] 
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The solution thus modeled is called Fowler and Guggenheim’s quasi- 
chemical model. 

This model may be calibrated on an experiment with a single adjustable 
parameter (3 q . Based on the knowledge of the values of the activity 
coefficients at a given temperature T, for a given composition, using 
relations [3.86] and [3.87], we calculate GT ; then, we use relation [3.85] to 
deduce J. 3 q and, with equation [3.79], we obtain the ratio Wp&lz - i.e. the 

difference £ m - £aa y £m . 

Let us calculate Warren and Cowley’s parameter for Fowler and 
Guggenheim’s quasi-chemical model. 

If we attempt to link our results to the values introduced in sections 3.1 
and 3.2.3, we immediately find: 


„\ _» _ 

z A — z B — z 

[3.88] 

(x) — ^aTaB — ^aTbA — ^B X A — ^A^B — ^AB ^ Z 

[3.89] 

^AB = ^BA =7 l=l ^ 

[3.90] 


A-*B 


By substituting this last relation into equations [3.76] and [3.81], we 
obtain: 


2zx a x b 
n = 1 ^-2- 

A-i 


[3.91] 


(3 q is defined, as usual, by relation [3.79], 

Thus, the activity coefficients expressed by relations [3.86] and [3.87] 
can be written in the form: 




1-x b (1-/7) 


[3.92] 
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]n/ B =^ln 


A 1-x a (1-/7) A 


[3.93] 


Figure 3.3 shows the variations of Warren and Cowley’s order parameter 
for Fowler and Guggenheim’s quasi-chemical model, for the following 
values of the parameters: T = 800 K, N a w AB = 30 kJ and z = 12. The result 
obtained respects the symmetry of the model with a minimum for the 
composition x A = x B = 0.5. 



Figure 3.3. Variation of the degree of order as a function of the composition of a 
binary solution in Fowler and Guggenheim ’s quasi-chemical model 
(reproduced from [DES 10], p.87) 


NOTE 3.8.- The values k AB and k BA defined by expressions [3.50a] and 
[3.50b] are li nk ed to one another by: 


(1-^7) = 


X A k AB 


X B k BA 


X B ( X B k AB X A ) X A ( X A k BA ~ X B ) 


[3.94] 


Thus, in Fowler and Guggenheim’s quasi-chemical model, the two 
weighting factors k AB and k BA depend on the composition. 

There are several models of solutions which use this concept of overall 
composition - Wilson’s, in particular. Given the restricted reach of these 
models, we have chosen to use the concept of local composition only in 
the context of the much more widely-applicable models UNIQUAC (see 
section 3.5) and UN1FAC (see section 3.6). 
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3.3. The quasi-chemical method of modeling solutions 

If we examine relation [3.31], it is exactly as though we had applied the 
law of mass action to a chemical equilibrium, which can be written as: 


A-A+B-B = 2 A-B [3R.1] 

The constant in this equilibrium would be equal to 1 , meaning that it is 
independent of the temperature. 

Similarly, by examining relation [3.75], we are led to write a law of mass 
action of the same equilibrium [3R. 1], but this time with an equilibrium 
constant of: 


K = exp- 


2w ab 
k J 


[3.95] 


It is because of these similarities that Fowler and Guggenheim’s model is 
known as the quasi-chemical model. 


Generalizations of this model have been put forward. Remembering that 
we have only taken account of the energies of interaction between near 
neighbors, we can extrapolate, by writing that the energy which is involved 
in the constant could be considered to be the sum of the energies between 
first neighbors, between second neighbors, etc., which would lead to a 
constant K appearing as the sum of the exponential terms, in the form: 


K = exp- 


2w ab 


+ exp- 


3w ab 


+ exp- 


4w ab 

vV 


+ ... 


[3.96] 


Thus, we can understand why the solution given by relation [3.31], which 
is the Bragg- Williams approximation, is also called the zero-order 
approximation, whereas Fowler and Guggenheim’s, given by relation [3.75], 
is called the l st -order approximation. Similarly, we can have 2 nd -order, 3 rd - 
order approximations, etc. 

As we can see, this quasi-chemical method is not dissimilar to the 
associated-solutions model discussed in section 2.5. This method has 
been rendered more generally applicable by writing what are known as 
quasi-equilibria - i.e. equilibria of the same type as above, but between 
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fictitious entities which might be atoms, molecules, vacancies, etc. - in fact, 
any entities the user cares to define. Thus, we can find pairs of molecules, 
trios, or even functional groups (parts of molecules), or clusters (groups of 
atoms or molecules). 

A generalization of Guggenheim’s quasi-chemical method has been 
constructed in order to find the number of possibilities of configurations 
between groups of elements for each pair in a binary solution of /V A 
molecules of component A and N B molecules of components B (here, the 
term “molecule” is used in a generic sense to speak of the elementary units 
involved in the processes). The units that come together in pairs may be true 
molecules, fragments of molecules, surfaces of molecules, parts of the 
surface of molecules, etc. The distinction is maintained between those 
elements which belong to component A and those belonging to 
component B. We suppose that the canonical partition function of mixing is 
of the form: 

zr=Zg(wA-)ex [3.97] 

k \ J 

g(a k ,b k ...) is the statistical weight of a configuration of elements in the 
volume defined by a series of particular values ct k ,b k of the configuration 
variables a, b, etc. These variables may be the volumic compositions, the 
local volumic compositions, the global or local surfacic compositions, 
depending on the elements in question. 

U mvc (a k ,b k ,...) is the configuration energy of mixing, which is a function 
of the configuration variables. 

Generally speaking, we do not know an exact solution to the problem. 
The quasi-chemical method gives an approximate value for our unknown: 
the total number of all the configurations. 

We know from Boltzmann’s law (see section 2.1.1) that in the sum 
involved in relation [3.97], one of the terms is much greater than the others, 
and that therefore the sum can be replaced by this maximum term. Hence, 
we can replace relation [3.97] with: 

[3.98] 


Z™ = g(a,b..)ex p 


U mix (a,b,...) 
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g{a,b..) and U ma ( a,b ,...) are the values of g and U mx which correspond to 
the most probable configuration - i.e. to the values a, b, ... of the 
configuration variables which yield the maximum Z"'“ . 

The hypothesis of the quasi-chemical method is to assume that the 
statistical weight g of the solution, i.e. the number of configurations with 
energy U mix that yields the maximum value of Z "' ,x , is given by the 
following equilibrium: 

Pair of elements belonging to A + Pair of elements belonging to B = Mixed pair 

The law of mass action applied to this quasi-chemical equilibrium is 
written as: 


g(a,b,...) = h(N A ,N B )g A (a,b,...)g B (a,b,...) [3.99] 

g A (a,b,...) denotes the number of configurations with energy U m,x 
which have the maximum number of homogeneous pairs A-A, and similarly 
g B (a,b,...) for the pairs B-B. The method postulates that h( N A , /V B , ) is 
contingent only upon the amounts of each of the components N\ and jV b 
present. The function h(N A ,N B ,) appears as a corrective factor, because the 
simple product g A g B would yield too high a number of configurations with 
the energy U. 

The solution must satisfy the boundary condition, which is that if the 
temperature increases indefinitely, U ma tends toward zero, and therefore we 
must find the number of configurations wherein the elements are distributed 
at random. For this purpose, we choose a reference model which satisfies 
that boundary condition. Guggenheim chose molecules comprising pairs of 
spheres of equal dimensions, distributed at random, which we examined in 
section 3.2.4, whereas in the UN1QUAC model, it is Staverman’s athermic 
solution for molecules of diverse shapes and sizes, distributed at random (see 
section 3.5) which is chosen. Suppose that we know the functions g A (a,b...) 
g B (a,b...) and g\a,b...) in this referential model. 
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To solve relation [3.99], we shall use a two-step process: firstly, we shall 
calculate the function h(N A ,N B ,), and then use the result in 
expression [3.99], 

In order to determine h{N A ,N B ,) , we must remember that this function 
depends not on the configurations but only on the amounts of the two 
components. Thus, we shall choose a solution of g that satisfies our 
boundary condition, which we know how to calculate. Clearly, the simplest 
such solution is that which we chose as a reference, for which the values are 
marked with a subscript a (e.g. g a ). This is the Bragg-Williams zero-order 
approximation. It obeys the following equation (random distribution): 

g a (a,b,...) = g aA (a, b,...)g aB (a, b,...) [3.100] 

As the chosen reference solutions are temperature-independent, the 
reference solution will be such that the function Z” IX in equation [3.100], 
when U"" x = 0 , is maximal. Thus, we shall write the system of equations 
which makes all the partial differentials of g in relation to the configuration 


variables equal to 0: 


d\ng a _dlng aA ( 31ng flB _ Q 

[3.101a] 

3 a k da k da k 

91ng fl _aing flA ( 31ng fl _ Q 

[3.101b] 

a b,. a b,. db. 


The resolution of this system leads us to the particular functions of the 
configuration variables a a (N A ,N B ) , b a (N A ,N B ) , ..., which enables us, on 
the basis of the values g A (a,b...), g* B (a,b...) and g* (a, b , ...), to calculate 
the particular functions g aA (N A ,N B ), g aB (N A ,N B ) and g a (N A ,N B ), and 
then, because h is independent of the model, we can use relation [3.99] to 
deduce the following from those functions: 

h(N A ,N B ) = 


[3.102] 
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Then, we move on to the second phase, searching for the so-called 
l st -order approximation or Guggenheim’s quasi-chemical approximation, 
characterized by the relation: 


f 

g (a,b,...)ex p 

V 


U' nix (a,b,...) 
R T 


g A (a,b,...)g B (a,b,...) 


[3.103] 


Now we shall cancel out the derivatives: 

ding U mix (a,b,...) _ d\ng A | dlng B _ Q 
da k R T da k da k 


[3.104a] 


ding U"" x (a,b,...) _ dlng A | dlng B 
d b k R T d b k d b k 


[3.104b] 


We deduce the functions a(N A ,N B ) and b(N A ,N B ), which enables us 
to calculate the functions g A (N A , N B ) and g B (N A , N B ) , and then by using 
expression [3.99], we can calculate the function g(N A ,N B ), which is the 
approximate solution to our problem. 

If the solution has more than two components, we write as many relations 
akin to [3.99] as there are pairs of components. For example, for a solution 
with three components, we shall have the following three equilibria: 

Pair A-A + Pair B-B = Pair A-B 

Pair A-A + Pair C-C = Pair A-C 

Pair B-B + Pair C-C = Pair B-C 

Hence, for A components, we would have /V(;V- 1 )/2 equilibrium states. 

Thus, the quasi-chemical method offers an approximate solution whose 
validity has been proven by Chang, with some extremely arduous work. 
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3.4. Difference of the molar volumes: the combination term 

In the models discussed hitherto, we have always assumed that the molar 
volumes of the two components of the solution were very similar, so that the 
imperfection of the solution was attributed solely to interactions between the 
molecules. Thi s hypothesis, though, is clearly not true in all cases - 
particularly in solutions of a polymer in a solvent - e.g. a solution of 
polystyrene in acetone, where the difference between the molar volumes is 
very significant: acetone has a molar volume of 73 cnrVmole, whilst that of 
polystyrene is, on average, 3333 cm 3 /mole. In order to take account of this 
fact, we introduce the concept of combinatorial excess entropy, which is 
linked to the distributions of the molecules in the space. 


3.4.1. Combinatorial excess entropy 

Flory and Huggins, simultaneously but independently of one another, 
constructed a term representing the excess entropy known as the 
conformation term or the combinatorial excess entropy. We shall now 
present their reasoning process for a mixture of small molecules, of a solvent 
A, and large one-dimensional molecules making up the solute B. The 
hypothesis at the heart of the model is still the pseudo-lattice, whose mesh is 
defined as follows: the molecule of component A (the smallest molecules) 
determines the acceptable dimension on each site of the lattice - its base 
volume. The larger molecules of polymer (component B) are virtually 
divided into sequences of the same volume as the small molecule, so that the 
large molecule contains V s sequences such that V. is equal to the ratio of the 
molar volumes of the two pure components: 


v = — 
' v,° 


[3.105] 


A molecule of B occupies the v s sites of the pseudo-lattice. If the 
solution contains Na molecules of component A (the solvent) and /V B 
molecules of component B (the polymer), the number of sites in the lattice 
is: 


N 0 = N a +v s N b 


[3.106] 



Microscopic Modeling of Liquid Molecular Solutions 93 


NOTE 3.9.- Nothing in the reasoning process suggests that the sequence of 
polymer chosen must correspond to a unit of monomer. The choice of 
sequencing relates only to the volume of the molecule of solvent. The 
sequence of polymer chosen may contain more or less than a monomer unit. 


We shall now construct the lattice in order to calculate the number of 
complexions £2, i.e. the number of possibilities for placing N\ molecules of 
A and N B molecules of B. Then, we shall use Boltzma nn ’s formula [2.41] to 
calculate the corresponding entropy term. 



Figure 3.4. Distribution of the molecules of solvent and polymer on the pseudo-lattice 


To construct the lattice in space, at each site, we place either a molecule 
of solvent or a sequence of a polymer chain, in the knowledge that those 
sequences are connected to one another. Figure 3.4 gives a 2D representation 
of the solution. Each grid square represents a site in the lattice which 
contains either a molecule of solvent A (in which case the square is white) or 
a sequence of the solute B (in which case the square is black). 

Each site is surrounded by z adjacent neighbors (z is the coordination 
number of the lattice). If we examine the closest neighboring sites of a 
sequence of polymer, there are two sites which must, inevitably, contain 
another sequence (if we discount the ends of the chains, of which there are 
very few, in relative terms). Hence, there are only z-2 neighboring sites that 
can be occupied by molecules of solvent A, and therefore able to interact in 
the mixture (interaction between sequences of polymer of two different 
chains also takes place in the pure polymer, and is not included in the mixing 
value). We suppose that the molecules of A and B are randomly distributed. 
Note that the entropic term of local composition which we encountered 
above is not brought into play; instead we restrict ourselves just to the 
combinatorial term. Thus, we have only one term of the excess Gibbs energy 
due to the distribution imposed by the volumes. 
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We shall break down the number of possibilities of placing /V B molecules 
of polymer on V S N B sites in a lattice containing N 0 sites. The remaining sites, 
obviously, are occupied by the Na molecules of solvent. 

Suppose that we have already put i molecules of polymer in place on the 
lattice. Let £X +I be the number of possibilities to place the (z+l) th molecule 
of polymer. We can calculate Q i+1 by considering each of its segments, 
one by one. Because the first i molecules are occupying iv s sites, there are 
N 0 —iv s places to put the first segment in place. The probability of a site 
closely neighboring this first segment being free is: 


Po = 


{N 0 -iv,-\) _ {N 0 -iv s ) 
K N 0 


[ 3 . 107 ] 


Because there are z sites neighboring this first segment, on average there 
are pz possibilities to place the second segment. That second segment has z- 1 
free neighboring sites, and the probability of those z-1 sites being vacant is: 


Pi = 


{K-K 

K 


2 ) 

-Po 


[ 3 . 108 ] 


Hence, the number of possibilities to place the third segment is p 0 (z - 1) . 
By the same reasoning, the same probability p 0 that the neighboring sites 
will be occupied is valid for the addition of any segment, provided we have: 

N 0 »v, [ 3 . 109 ] 

Thus, the number of possibilities of placing the v s segments of the (z+l) th 
polymer molecule on the lattice is: 

n i+l ={N 0 -iv s )p v 0 ‘- l z(z-l) v ’- 2 [3.110] 


If z is not too small and V s is sufficiently large, we can replace 
z(z - l) v '“ 2 with (z - 1) 1 '* -1 . By doing so, we obtain: 


^,+1 = { N 0~ iV s ) Po~' z ( z - 1 )"' s “ 1 


[ 3 . 111 ] 
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In light of expression [3. 107] for p 0 , this gives us: 




N n 




noting that with the condition [3.109], we can write: 


[3.112] 


[7V 0 -(/ + lv s )]! [ 3 .H 3 ] 

(N 0 - iv s )(N 0 - iv s - 1)...(N 0 -(/ + \)v s + 1) = (N 0 - iv s f 


By substituting this back into equation [3.1 12], we find: 

[3.114] 


[3.115] 

However, the number of complexions or number of possibilities of 
placing the N B molecules of polymer will be such that: 

Q = TnIl Q m P-H6] 

■'’B ' 1=1 

The division by jV b ! stems from the indiscemibility of the molecules of 
polymer. 

If we feed expression [3.115] back into relation [3.116], we can separate 
the product into two groups in the form: 


r> ~ (*o *,)'■ z l.^-i 

,+1 = [1V 0 -(i + 1)v,]! C 1V 0 ) 
This can also be written as: 

[W,-(,-l)K,]! z -i 
(W 0 -,V,)! 


1 

f Z_1 ] 

JVb(v.-I) 

A^o-O'-ik]! - 

N ! 

iV B * 

V ^0 J 


1 

T 

£ 

H ~ 

1 


Q 


[3.117] 
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We can show that the second factor in the above relation can be 
written as: 


' ■ b 

n 


|X-0'-lk]!_ N 0 l 


(^o"*V,)! {N 0 -N b v,)\ 


[3.118] 


Thus, the number of complexions will be: 

Wb(v s -1) 


n = 


(N a +N b v s )\ 


f 


7V a !AV 


Z — 1 


\N A + N B v sJ 


[3.119] 


By switching to logarithms and applying Stirling’s approximation, we 
find: 


. „ A , , N a + N b v N a +N b v s 

InQ = N a In — — + N 7 In — — 

N, N 2 

+tf B (v,-l)[ln(z-l)-l] 


[3.120] 


Thus, if we substitute equation [3.120] back into Boltzmann’s relation, 
we find the entropy term of conformation: 


S (C) = R 


- N B 

+«b - !) [ ln (^ - !) - 1] 


[3.121] 


If we denote the volumic fraction of component i as: 
N, 


0 = 


N a +N b v s 


[3.122a] 


the conformation entropy can be written as: 

S lC) =- R {« A ln^ A +n B \n0 B ~n B (v s -l)[ln(z-l)-l]} 


[3.122b] 
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However, the conformation entropy of the pure solvent on the lattice is 
zero (there is only one configuration with one molecule per site), and that of 
the polymer is given by relation [3.121] with Na = 0 - i.e.: 

S p = R« b {In v, + (v s - 1) [ln(z - 1) - 1]} [3. 123] 

Thus, the conformation entropy of mixing of the solution is: 

S mix(C) = 5 (C) -S p = -R{h a ln<Z> A + n B ln<Z> B } [3.124] 

If we subtract the entropy of mixing from the perfect solution given by 
relation [A.2.10], we obtain the following for the molar excess conformation 
entropy: 


S; (C) =-R r A ln^+x B ln^ 

x R 


[3.125] 


We use this expression for a model of an athermic solution, defined in 
Chapter 2 (see section 2.6). 


3.4.2. Flory’s athermic solution model 

Flory’s model applies to an athermic solution, meaning that the excess 
molar enthalpy is null, so the excess molar Gibbs energy is simply written as 
follows, in light of relation [3.125]: 

G: =~TS: =R T {x A ln^A + x B ln^4 [3.126] 

I *A J 


For a binary solution, we use the excess partial molar Gibbs energy to 
deduce the activity coefficient of a component i: 


m 0 0 V V 

In r ] 7 = In — L + 1 — ^ = ln^ + l — i- 

X,. x,. V„, F 


[3.127] 
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V, 

Tl (/) 

Tb (/> 

10 

0.245 

0.0012 

100 

2.7xl0‘ 2 

10‘ 4U 

1000 

2.7xl0' 3 

o 

i- 

o 


Table 3.1. Values of activity coefficients with infinite dilution of a polymer 
(B) in a solvent (A) as a function of the “length ” v s of that polymer 


In particular, at infinite dilution, we obtain the activity coefficient of 
each of the two components of a binary by making n B = 0 for that of 
component A and n A = 0 for that of component B. This gives us: 

In = ln-k + 1 — -f = In— + 1 — — [3.128] 

v 2 v 2 ° v, v. 

In r ; (I) = In 4 + 1 - 4 = K + 1 " K [3-129] 

v, v. 

In this model, we have not taken account of the energetic interactions 
between the molecules. For precisely this reason, such as model is well 
suited for solutions of polymers, as the intermolecular energy values are very 
low indeed in these solutions. Table 3.1 shows that the deviation from ideal 
state grows very rapidly as the “polymer chain length” v s increases. 


3.4.3. Staverman ’s corrective factor 

In the above model, we chose linear polymer molecules and molecules of 
solvent which are much smaller than those of the polymer, so that every part 
of the polymer was in the vicinity of molecules of solvent. 



Figure 3.5. Exclusion of sites available for the solvent 
due to the closure of the polymer molecule 
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Staverman examined the problem of representing a solution in which the 
molecules are no longer linear, and therefore a certain amount of the 
molecule would not be exposed to the molecules of solvent, because they 
would be isolated inside the molecule of solute. These imprisoned sites 
receive no molecules of solvent. Figure 3.5 gives a 2D representation of the 
exclusion of certain sites (shaded) which are no longer able to receive 
molecules of solvent as near neighbors. Thi s introduces a corrective factor 
relating to the surface of the molecule. 

3.4.3. 1 . The concept of structural parameters 

We shall now introduce two fundamental parameters of the molecules 
which are used in the models of molecular solutions discussed below. These 
two parameters are: 

-the ratio of the van der Waals volume of the molecule to the van der 
Waals volume of a standard molecule or segment: 


v 

_ v «■(/) 

v 

v w(st) 


[3.130] 


- the ratio of the area of the molecule to the area of a standard molecule 
or segment: 


9 t = 


■^w(st) 


[3.131] 


The choice of standard is arbitrary. Frequently, we use the Abrams and 
Prausnitz standard [ABR95]: a linear molecule of polyethylene of infinite 
length, which conforms to the following equation: 


T f - 1 = |( r „ where z = 10 [3.132] 

This gives us: 

V w(sl) = 1 5. 1 7 cm 1 /mole and H w(st) = 2.5x1 0 9 cm 2 /mole 


[3.133] 
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Thus, the structural parameters are calculated on the basis of the two 
relations: 

r, =V w(i) / 15,17 and q,=Z Mi) ! 2.5xl0 9 [3.134] 


Table 3.2 gives the values of the structural parameters of a several 
molecules, calculated by Abrams and Prausnitz on the basis of the van der 
Waals volumes and the areas published by Bondi in 1968. 


Based on the structural parameters, we define: 

- the volumic fraction of a component i by: 

. Nr xr 
0 = — — = — — 


a I 


xr 


- the areic fraction of a component i by: 

N i<li x t q, 


e,= 


Z N i c ti Z v 3/. 

i = 1 1=1 


[3.135] 


[3.136] 


Fluid 

r i 

9; 

Water 

0.92 

1.40 

Carbon dioxide 

1.30 

1.12 

Ethane 

1.80 

1.70 

Benzene 

3.19 

2.40 

Toluene 

3.87 

2.93 

Aniline 

3.72 

2.83 

n-Octane 

5.84 

4.93 

Acetone 

2.57 

2 .34 

Dimethylamine 

2.33 

2.09 

Acetaldehyde 

1.90 

1.80 

Polyethylene 

0.67v, 

0.54v, 


Table 3.2. Values of the structural parameters for various molecules 
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3.4.3. 2. Staverman ’s model 

Staverman uses Flory’s model as a starting point, but introduces a form 
factor in the number of complexions, taking account of the fact that the 
molecules are only in contact at their surfaces. This enables him to construct 
a new expression for number of possibilities of introducing the molecule of 
polymer and the molecules of solvent. Thus, the excess molar entropy term 
of conformation is altered, and becomes: 




[3.137] 


We can see that the form factor is illustrated by the ratio 0. / . 

Using the concepts which we have introduced, we shall now present one 
of the most popular models of solution: the UN1QUAC model. 


3.5. Combination of the different concepts: the UNIQUAC model 

The UNIQUAC (Universal Quasi-Chemical) model was introduced by 
Abrams and Prausnitz (1975), using Guggenheim’s quasi-chemical model 
and applying the concepts of conformation with Staverman’s relation and 
Wilson’s local-composition model. 

This model, which yields excellent results for polar and non-polar 
molecular liquids, is especially well suited for the study of liquid/ 
vapor equilibrium and the equilibrium between two liquids that are not 
completely miscible. Regardless of the number of components of the 
solution, the application of this model only requires the knowledge of two 
adjustment parameters per binary system, which can be deduced from the 
solution. The model is so widely applicable that it actually contains a 
number of previously classic models such as the models put forward by 
VanLaar, Wilson, Renon et al. (the NRTL - Non Random Two Liquids - 
model), Scatchard and Hildebrand, Flory and Huggins as special cases. 
In addition, it lends a physical meaning to the first three coefficients 
/3 j ,S j and £ t in the Margules expansion (equation [2.1]). 
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The model is founded on the following hypotheses: 

- we accept the existence of a pseudo-lattice in space, which occupies the 
whole of the volume. The vacancies in the liquid will not be taken into 
consideration, either in the pure liquids or in the solution; 

- the coordination number of the lattice is the same around an element of 
A as it is around an element of B, and it is also the same as in each of the 
pure liquids of which the solution is composed. Thus, for example, for a 
binary solution, we would have: 

z a= z b= z \= z 'b [ 3 - 138 ] 

-the interacting elements chosen are segments of the surface of each 
molecule. Each molecule may present different surface elements; 

- the fragments of components are chosen in such a way that their 
dimensions, volume and surface area are practically the same. As we saw 
above that vacancies are not taken into account, this means that the mixing 
volume is null, and therefore that we can treat the internal energy and 
enthalpy of mixing as one and the same thing, and do likewise for the 
Elelmholtz energy and the Gibbs energy of mixing; 

- although the interacting surface elements may be different in nature 
from one point in the solution to another, we shall make the hypothesis that 
the energies of interaction are mean values, which depend only on the origin 
of the elements. Thus, for instance, for a binary system, we would have the 
energy £ BA = £ AB between two surface elements, one of which belongs to the 
molecule of component A and the other to the molecule of component B, 
£ aa between two elements belonging to two molecules of component A and 
£ bb between two elements belonging to two molecules of component B. 

The states will be characterized by the local areic fractions: 

d AA denotes the surface fraction of surface elements derived from A 
around a surface element derived from A. We can immediately deduce the 
mean of the other symbols: 0 m , 0 AH and 0 BA . Of course, the balance 
around any surface element gives us: 



Microscopic Modeling of Liquid Molecular Solutions 103 


^BB ^AB 1 


[3.140] 


Thus, as configuration variables, we can choose a value belonging to each 
of the two pairs 0 AA or 0 BA , on the one hand, and 0 BB or 0 AVi on the other. 

In light of our last hypothesis concerning the energies of interaction, the 
internal energy of mixing will be given by the expression: 


U ma - u;r = + 0 BA £ BA ) 

_ (^BB^BB ^AB^AB ) 


For the sake of ease, we shall set: 


[3.141] 


11 a = 2 e a ~ u a [3-142] 

This gives us the following expression for the internal energy of mixing: 

[3.143] 


U ma ~U 0 = -q A N A {0 AA u AA + 0 BA u BA ) 

_ ?B^B (^BB M BB T ^AB M AB ) 


We shall express the canonical partition function of mixing by using 
Guggenheim’s quasi-chemical method (see section 3.3), applied with 
Staverman’s athermic model as a reference, where the excess entropy is 
given by relation [3.137]. 

In Staverman’s model, the statistical weights of the homogeneous 
elements g* and g* are given by: 


8i = 


§ 2 = 


{q l N i 8 n + q 2 N 2 6 n ) ! 

{q^AMq^Y- 

(q 2 N 2 0 22 + q l N x G n ) ! 
(q 2 N 2 6 22 ) \(q 2 N 2 O n ) ! 


[3.144] 


[3.145] 
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With the zero-order approximation, calculation leads us to: 


>3 

II 

> 

II 

S K 

qb 

II 

[3.146] 

# A-^A 



n(0) ‘Zb^B 

t'R A _ 

_ n( 0) _ a 
— &BB — ^B 

[3.147] 

QaN A T ‘Zb^B 




0 j is the areic fraction of the molecule i. It is given by relation [3.136]. 
The application of the l st -order approximation yields: 


6aa =- 


0 A 


0 A + 0 H eX P 


f u ba~ u aa^ 


R T 


[3.148] 


^bb =- 


0» 


0 B + 0 A eX P 


^ U AB U BB ^ 

R T 


[3.149] 


After calculating the partition function using relation [3.97], we deduce 
the excess molar Gibbs energy using the relation: 


XS xs 


RTTnZ 


[3.150] 


We can see that this molar Gibbs energy can be expressed as the sum of 
two terms: 


Q xs ^xs(Com) 


-rX5(Res) 


[3.151] 


One of them, 1 9 i s drawn from Staverman’s athermic solution, 

from which it draws its name - the conformation term : 


'fxs(Com) 


= R T 


0 0 

x. In — — + x a In — - 


z 

+ — 
2 


, 0 A 0R 

<lA x A^-ff+q B x^n-f- 


V 


<p. 


0 , 


B 7 


[3.152] 
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The other term, G„“ (Rcs) , is called the residual term. Its value is: 

G,f ResJ = RT[-g A x, \n(e, + e H T HA )-q R x H ln(0 B +0,r AB )] [3.153] 

In this latter expression, the different symbols have the following 
meanings: 

[3.154a] 





(>; - 4 , )-(>;- !) 


[3.157] 


with a symmetrical expression for component B. 
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In the view of the model’s developers, any value of the coordination 
number z between 6 and 12 has practically no impact on the results. Thus, 
they chose the value: 


z = 10 


[3.158] 


For a solution with N components, the above results can be extrapolated, 
and we obtain the following results: 


■^xs(Com) 


= R T 


N _ ( N A 

Z*iln — + T 


V ,=1 


0, 


ij 


[3.159] 


■'rxs'(Res) 


R T 




7=1 


T n = exp 


f Uj t - uii ^ 
R T 


= exp 


4 a .. ^ 

j' 

V T 


[3.160] 

[3.161] 


The activity coefficients of each of the components are also divided into a 
conformation term and a residual term, and we have: 


]n V / )(Com) =]n ^ + £ ? ]n i + / 

X; 2 ‘ 0 ' X ^ , J 


7=1 


[3.162] 


In = q. In 


N N a T 


k^kj 

V t=l J 

M) _ v (7)(Com) , , ,//)( Res) 


In /. 1 = In y) } > + In y.‘ 


[3.163] 


[3.164] 


We can see that the UN1QUAC model requires only two adjustable 
parameters for each pair of components: the differences: u tJ - u /; and w ,, — u a 

(or dj and a ..). However, we shall take the precaution of not involving too 
many parameters. For a system with N components, there are only 
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N(N — l)/ 2 binary systems and therefore the same number of pairs of 
adjustable parameters. 


Name of model 

Approximations about 
equation [3.152] 

Approximations about 
equation [3.153] 

Flory-Huggins 
athermic model 

q x - r\ and q 2 = r 2 

U 2| — ^11 — ^22 

Wilson 2-parameter 
model 

q x = r x and q 2 - r 2 

ii 

CN 

II 

c? 

NRTL model 

✓~rXS(Com) rv 

KJ m 

(^l) (^ 2 ) ^12 

Van Laar’s model 

xs(Com ) rv 

v ~m 

Expansion of (zk/ / RT) , limited 
to the 3 rd term 

Scatchard- 
Hildebrand model 

✓~rxs(Com) rv 

m 

Expansion of (zlw / R T ) limited 
to the 3 rd term 

Margules 3-tenn 
model 

s~ixs( Com) /v 

m 

Expansion of (zkt / R T ) limited 
to the 4 th tenn 

Flory-Huggins non- 
athermic model 

q x = r\ and q^ - r 2 

Expansion of (Aw / R T ) limited 
to the linear tenn 


Table 3.3. Particular models of solutions included in the UNIQUAC model 


The model’s creators showed that a certain number of classic models of 
solutions published previously were actually included in the UNIQUAC 
model, as they appear to be particular cases. By way of certain 
approximations performed on equations [3.152], [3.153] and [3.154], we 
obtain the expressions of the excess molar Gibbs energy for these different 
models. Table 3.3 lists some of these simplifications. 


3.6. The concept of contribution of groups: the UNIFAC model 

The UNIQUAC method requires complete calculation in each specific 
case. Now, though, we are going to examine the concept of the contribution 
of groups, which can be used to calculate certain properties of the molecules 
a priori, using relatively small databases, and the application to a model to 
calculate the excess Gibbs energy a priori. 




108 Modeling of Liquid Phases 


3.6.1. The concept of the contribution of groups 

The concept of the contribution of groups is based on the idea of 
additivity. Some of the properties of a molecule can be expressed by a 
weighted linear combination of the contributions of the groups of atoms 
making up the molecule. 

Thus, for example, we can calculate the values of the parameters r. and 
q i of the volumic fraction and areic fraction of a molecule, as being the 
weighted combined of values of those same properties attributed to groups of 
atoms. Thus, for the volumic fraction parameter r. , we would write: 

r, [3-165] 

k 

vf is the number of groups of nature k present in the molecule i, and R k 
is the parameter of the volumic fraction attributed to group k (defined on the 
basis of equation [3.134], as was the case for molecules). The sum is 
extended to all the groups of which the molecule is composed. 

Similarly, for the parameter of the areic fraction q t , we shall have: 

[3.166a] 


To do this, we have divided the molecules into groups of atoms, which 
are, themselves, divided into subgroups. Each subgroup is assigned a value 
R k and a value Q k . Table 3.4 gives the contributions of a number of groups. 
Each first term of an organic series constitutes a group. 

The values calculated by relations [3.165] and [3.166] can be used in the 
UN1QUAC model (see section 3.5) if the molecular values are not available. 


3.6.2. The UNIT AC model 

Models known as group models have been created, based on the 
following observation: the number of systems - even just the number of 
binary systems - is extremely high, and new ones are being discovered all 
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the time; meanwhile, the number of functional groups of which these 
molecules are composed is far less. 


Group 

Subgroup 

A 

Qk 

Cl 1 

CH, 

0.9011 

0.8480 

CH : 

0.6744 

0.5400 

CH 

0.4469 

0.2281 

C 

0.2195 

0 

c=c 

CH 2 =CH 

1.3454 

1.1760 

CH=CH 

1.1167 

0.8670 

ch 2 =c 

1.1173 

0.9880 

CH=C 

0.8886 

0.6760 

C=C 

0.6605 

0.4850 

ACH 

ACH 

0.5313 

0.4000 

AC 

0.3652 

0.1200 

OH 

OH 

1.0000 

1.2000 

CH,OH 

CH,OH 

1.4311 

1.4320 

ACOH 

ACOH 

0.8952 

0.6800 

CH 2 CO 

CH 3 CO 

1.6724 

1.4880 

CH 2 CO 

1.4457 

1.1800 

CHO 

CHO 

0.9980 

0.9480 

COOH 

COOH 

1.013 

1.2240 


Table 3.4. Functional groups and subgroups and structural parameters for the UNIFAC 
model (“AC” represents a carbon belonging to an aromatic ring) 

The UNIFAC model (UNIQUAC Functional Group Activity 
Coefficient), which was put forward by Fredenslund, Jones and Prausnitz in 
1975, is a group model where the idea is to use the existing data 
on equilibrium states to predict the properties of systems for which 
no experimental data are available. Thus, it is an entirely predictive model 
which, unlike the models discussed above, does not require us to determine 
parameters by calibrating the experimental results found for partial systems 
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(in contrast to the UNIQUAC model, which requires the data 
from binary systems in order to deduce the properties of more complex 
systems). 

The group-model method requires us to work back from experimental 
results to determine parameters which are characteristic of the interactions 
between pairs of structural groups, and then inject that knowledge back into 
models to obtain the properties of new systems containing those same 
structural groups. 

The idea underlying the UNIFAC model is the consideration that a 
solution of molecules i, j, etc. behaves like a solution of the functional 
groups k, m, etc. making up those molecules. 

For example, a mixture of linear alkane molecules involves three types of 
interactions: CH 3 -CH 3 , CFfr-CFF and CH 2 -CH 2 . The properties of any given 
mixture of alkanes must be able to be deduced from the properties linked to 
the three types of interactions at play. Of course, this whole construct is 
founded on a hypothesis which stipulates that the interactions between 
groups do not depend on the environment of these groups in their respective 
molecules. Yet, indubitably, it is this hypothesis which is the main source of 
errors that become apparent when the calculation results are contrasted with 
experimental data. For example, we can see that, in a complex system, the 
results given by the wholesale use of UNIFAC are much poorer than those 
given by the use of the UNIQUAC model. For this reason, the UNIFAC 
model is not used wholesale to model a complex system; instead, essentially, 
it is used to obtain the necessary data on binary systems for which no 
experimental data are available. These data are then fed back into a 
UNIQUAC-type model (see section 3.5). 

Thus, we shall apply the UNIQUAC model to this solution of functional 
groups. 

From this, we deduce that the activity coefficient for a component i obeys 
relation [3.164], meaning that its logarithm will be the sum of two 
contributions - one conformation contribution, pertaining to the entropy, due 
to the arrangement of the molecules within the volume, and the other 
residual, relating to the enthalpy and due to the different interactions 
between pairs of molecules. 

The conformation contribution is written by following relation [3.162] - 
indeed, the arrangement of functional groups obeys the same laws as that of 
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entire molecules, in view of the unavoidable intramolecular bonds between 
the functional groups in the same molecule. 

The residual contribution is obtained by considering a functional-group 
solution. 

The values corresponding to the functional groups will be written with 
the same letter as the same value in a solution of molecules, but in 
uppercase. Thus, the coefficient of activity of the group k will be r k , its 
volume parameter R k , its surface parameter Q k , its molar fraction X k . In 
the energy exponential of relation [3.153], the letter T would be replaced by 
H* km , its surface fraction 6 would be replaced by 0 k . Only the term a in 
the equivalent relations [3.154] would retain the same notation, a hn . By 
transposition of relation [3.163] to the group solution, the residual activity 
coefficient r k for a group k will be written as: 


RT\nr k =Q k 



z 


Z 6 ^ 


[3.166b] 


The different entities are defined, as before, for the groups this time by: 
The surface fraction of the group k\ 


& k = 


x k Q k 


Z x k Q k 

k 

The energy term between the groups k and m: 

'f' km = eX P f 


[3.167] 


( a k,n ^ 

T 


[3.168] 


The molar fraction of the group k\ 

IT- 


X, 


x k = 


zz 


v k ) x i 


[3.169] 
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Of course, that same coefficient of residual activity of the group is, by 
definition: 


RT In r t = n k - nl [3.170] 

In this expression, /u k is the chemical potential of the group k in the 
group solution studied, and ju k is the chemical potential of component k in a 
hypothetical substance in which the component k is pure. Evidently, such a 
substance can neither be formulated nor conceived. However, we are able to 
conceive a pure solution of i, considered to be a solution of its groups, in 
which the group k has a chemical potential jj { k ] and we can write the 
following for the chemical potentials and the residual activity coefficients: 

<u k -rfP=(M k ~K)~ (/4° -K) = Rr(lnr i -lnrf ) [3.171] 

r[:' is the coefficient of residual activity of the component k in a pure 
solution of “component” i. 

In addition, the residual chemical potential of component i is given by the 
linear combination of the chemical potentials of the different groups, in 
accordance with: 

/A =IXX [3.172] 

k 

Similarly, for the pure solution of i: 

[3-173] 

k 

This gives us the coefficient of residual activity of component i (using 
convention (I) - pure substance) in the solution at hand: 

RTTn y[ Res) = A - g,° = RT^ (in - In ) [3.174] 

k 

The term In r k is calculated on the basis of relation [3.166]; the term 
In/’) 0 is also calculated on the basis of the same relation [3.166], but 
introducing the elements for pure i - in particular, the molar fraction of the 
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group k is in the pure component i and is given by expression [3.169], 
making x, = 1. This term is then fed back into relation [3.167] and the 
interactions to be taken into account in expression [3.166] exist only 
between groups belonging to the component i. 


Groups 

ch 2 

ACH 

ACCH, 

OH 

CH 3 OH 

ACOH 

CH 2 CO 

CHO 

COOH 

ch 2 

0 

61.13 

76.5 

986.5 

697.2 

0133.3 

476.4 

677 

663.50 

ACH 

-11.12 

0 

167 

636.1 

637.3 

1329 

25.77 

347.3 

537.4 

ACCHt 

-69.7 

-146 

0 

803.2 

603.3 

884.9 

-52.1 

586.8 

872.3 

OH 

156.4 

89.6 

25.82 

0 

-137.1 

-259.7 

84 

-203.6 

199 

CHjOH 

16.51 

-50 

-44.5 

249.1 

0 

-101.7 

23.39 

306.4 

-202.0 

ACOH 

275.8 

25.34 

244.2 

-451.6 

-265.2 

0 

-356.1 

-271.1 

408.9 

CH 2 CO 

26.76 

140.1 

365.8 

164.5 

108.7 

-133.1 

0 

-37.36 

669.4 

CHO 

505 .7 

23.39 

106.0 

529 

-340.2 

-155.6 

128 

0 

497.5 

COOH 

315.3 

62.35 

89.86 

-151 

339.80 

- 11.00 

-297.8 

-165.5 

0 


Table 3.5. Energy interaction terms of groups, expressed in Kelvin 1 
(“AC” represents an aromatic carbon) 


All of the activity coefficient of the component i in the solution, 
therefore, is given by relation [3.164], in which the terms are given by 
relations [3.162] (with [3.156]) and [3.173] (with [3.166]). 


We can easily calculate the excess molar Gibbs energy and the excess 
molar enthalpy, using: 


G, 


xs 

m 




[3.175] 


H' 


=-rt\YYxv 


k 


T 


3'nr, 
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-T 


J P,n 


ainr'"') 


dT 


J P,n 


[3.176] 


The calculations are conducted using values of the structural parameters 
of the groups and the energy terms of pairs of groups a km . The initial data, 
drawn from the study of numerous liquid/vapor and liquid/liquid equilibria, 
were provided by Gmehling in 1979. A few examples of structural 
parameters and energy terms a km (in Kelvin -1 ) between two groups k and m 
are listed in Tables 3.4 and 3.5 respectively. 
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3.6.3. The modified UNIFAC model (Dortmund) 


The UNIFAC model discussed in the previous section does suffer from a 
number of shortcomings, particularly for values of the activity coefficients in 
infinite dilution, the enthalpies of mixing and the values obtained for highly- 
asymmetrical systems. In order to improve the data, adjustments were made 
in three directions (1987): 

- a refinement of the functional groups, and therefore of their number. 
The UNIFAC Consortium (Dortmund Data Bank), a user group, now 
produces 92 groups instead of the 67 that were included in the first system. 
The result of this is that we are better able to take account of the groups’ 
environments; 

- an increase in the number of binary systems studied experimentally, 
with a view to refining the energy values; 

-better integration of the variation of the results with temperature, by 
writing the energy parameter in the form: 


A J- 


' C kn,T 


[3.177] 


Obviously, the consequence of this is the provision of four extra 
parameters per pair ( b km , c km , b mk and c mk ). 


All of these improvements were brought together in the creation of a new 
model, called the Dortmund modified UNIFAC model, which is 
characterized, for each binary system, by four structural parameters and six 
energy parameters contained in a database. 

The UNIFAC consortium also gives its subscribers a piece of software 
capable of calculating various thermodynamic values on the basis of the 
model and its database. 


3.6.4. Use of the UNIFAC system in the UNIQUAC model 

Although the UNIFAC model is totally predictive, it is rarely used for 
systems with more than two components because, as indicated above, the 
results given in such cases by the UNIQUAC model are often of better 
quality. On the other hand, UNIQUAC requires us to know two parameters 
per pair of components: the energy parameters, drawn from experimentation. 



Microscopic Modeling of Liquid Molecular Solutions 1 1 5 


Occasionally, when using it, we do not have the necessary data for all the 
pairs that it is necessary to consider. Therefore, we use the UNIFAC model 
to provide the data corresponding to certain pairs for which we have no 
experimental data. It is then necessary to transfer parameters from the 
UNIFAC model to the UNIQUAC model. This is what we shall now 
examine. 

Remember that the use of UNIQUAC requires parameters T km and 
which are at play in relation [3.153] - in fact, by relations [3.154a] and 
[3.154b], the energies a km and a mk . It should be noted that, in spite of the 
similitude of the notations, these values are not those that are given by the 
group data (e.g. in Table 3.5) which are group values, whereas those 
required by UNIQUAC are mean values per molecule. In order to switch 
from one to another, we proceed as follows. 

For the binary system in question, we calculate the values of R^lny; tRes, 
for each of the two components, on the basis of the UNIFAC system and 
relation [3.174], for different values 0 n i.e. the composition of the binary 
solution defined by the quantities NA and NB of molecules of the two 
components and relation [3.136]. We then adjust the two curves 
RTTn j4 Res >(6/ A ) and R7Tn ;4 Res> (0 B ) on the UNIQUAC model on the 

basis of relation [3.161] by choosing appropriate values Tab and x BA . These 
latter values, or the corresponding energies a AB and a BA , are then used with 
UNIQUAC to study the complex system containing the pair A-B. 

In this chapter, we have examined a number of models based on what is 
known as the G xs method. They share a common characteristic: they do not 
take account of any possible vacancies in the liquid, i.e. non-null excess 
volumes. Another method, based on equations of state, gives better results in 
certain cases, particularly when one of the components of the solution is in a 
hypercritical state and in the vicinity of the critical region of the solution. 
These models are described in the chapters devoted to mixtures of gases. 




4 


Ionic Solutions 


Ionic solutions are set apart from molecular solutions by the intervention 
of Coulombian forces between electrical charges - forces which decrease far 
more slowly with increasing distance than do the van der Waals forces 
which exist between molecules. An important consequence of this difference 
is that for ionic solutions, we can no longer content ourselves with the two- 
body interactions between near neighbors. This means that, in the case of 
ionic solutions, the approximation of the ideal dilute solution model 
(see sections 2.4 and 3.1.3) is acceptable only for much lower concentrations 
than in the case of molecular solutions. In other words, for an identical level 
of concentration, an ionic solution is more imperfect than a molecular 
solution. 

Numerous models of ionic solutions have been put forward in the 
existing body of literature. The most important of these models, which is 
actually found to be included in all the others, is Debye and Hiickel’s, 
which attributes the imperfection solely to the electrostatic forces 
between the ions but, in spite of this, is acceptable only for fully-dissociated 
strong electrolytes and very-dilute solutions. Then, we shall cite 
Pitzer’s model (1973), which combines Debye and HiickeTs model with a 
virial-type expansion, and is therefore able to extend the range of 
concentrations examined. Beuner and Renon’s model (1978) builds on 
Pitzer’s, extending it to solutions containing neutral molecules such as SO 2 , 
NH 3 , CO 2 or H 2 S. The most recent models take account of the concept of 
local composition (see section 3.2). In this category, we can cite the NRTL 
electrolyte model introduced by Chen (1979). Finally, other models have 
expressed the intermolecular interactions over short distances by a 
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UNIQUAC simulation (see section 3.5), and even the UNIFAC version (see 
section 3.6). 

The objective of all the improvements made to modeling has always been 
to find models that cover ranges of increasingly high concentrations, which 
are usable both for weak electrolytes and strong electrolytes, i.e. introducing 
models of solutions that are mixtures of ions and molecules. 

In this context, we can agree that the interaction energy can be considered 
to be the sum of two contributions, in the form: 

w = mA'' + w (M oI) [4.1] 


The first term w (d> is a contribution due to the electrical charges, 
connected to long-distance interactions; the second term if Mo1 ' corresponds 
to intermolecular forces over a short distance and pertains to the molecules 
and the ions independently of their charges. 

In general, the first term yields Debye and Hiickel’s model (see 
section 4.2). 

The second term can be expressed in a variety of forms, such as the form 
of a virial involving the two-way and three-way interaction terms of the 
molecules (etc.) as found in Pitzer model (see section 4.3), and yield an 
excess Gibbs energy in the form: 

Q xs ^rxs(elec) _j_ ^xs(Virial) |-^ 

Another possibility is that the excess Gibbs energy can take the form of a 
sum of two terms such that: 


~ixs /'■fxs’(elec) 


+ G 


«(UNI) 


[4.3] 


- a term of electrostatic interaction, the basis for which is still Debye and 
HuckeTs model, which pertains only to the ions; 

- a “UNIQUAC” or “UNIFAC” term, concerning both the ions and the 
molecules, which is of the same nature as that which we encountered in the 
molecular models (see sections 3.5 and 3.6); 
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Relations [4.2] and [4.3] are transposed to give relations [4.4] and [4.5] in 
terms of activity coefficients: 


In ft = In yf elec) + In yf VMal) 

[4.4] 

In y i = In yf elec) + In yf UNI) 

[4.5] 


In what follows, we have chosen to show the progression of the modeling 
of ionic solutions, focusing at length on Debye and Hiickel’s essential 
model, Pitzer’s model (which is very widely used in industrial circles) and 
the extension of the UNIQUAC and UNIFAC models to mixture of ions and 
molecules. 


4.1. Reference state, unit of composition and activity coefficients of ionic 
solutions 

The molecular models that we have discussed (see Chapters 2 and 3) all 
yielded activity coefficients in reference (I), pure substance reference, and 
the compositions were always expressed in molar fractions. In the case of 
ionic solutions, firstly the distinction between solvent and solute is often 
very clear, and secondly as the notion of a pure ion is illusory, usually the 
activity coefficients of the components of an ionic solution are expressed in 
reference (II), the infmitely-dilute solution, or even more frequently, in 
reference (III), the molar concentration. As is suggested by relations [4.4] 
and [4.5], we shall add terms put in place in the case of molecular solutions 
and purely ionic solutions. Thus, in order to perform this addition, we need 
to be able to switch, for the molecular components, from an activity 
coefficient in reference (I) to an activity coefficient in reference (II) or (III). 
Relation [A.2.21] is used to switch from jA°to , and given that we 
can write: 


Ms ~g, 

R T 


= ln x! 




[4.6] 


relation [A.2.21] is transformed into: 

rl U) = y? I /r 


[4.7] 
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This relation will enable us to add terms of the same nature into a sum of 
activity coefficients such as relation [4.4]. As we can see, the electrostatic 
term is given directly in convention (II), whereas the combination and 
residual terms are, as we have seen (in section 3.5), obtained in 
convention (I). 

Additionally, in normal usage, the composition of an ionic solution is 
very frequently expressed in terms of molality ( M s : number of moles per kg 
of solvent) or concentration (moles per liter of solution), rather than molar 
fraction. 

We find the molar fractions on the basis of the molalities using the 
relation: 


k A ^ K r A OT 

M = — = — — = — — [4.8] 

m 0 n 0 M 0 x 0 M 0 

The chemical potential of the component 5 becomes: 

A = M: + R T In / S U) X 0 M 0 M S = //( + R T In M 0 + R T In M s / s m) [4.9] 

We find the following equivalences, in light of relation [4.8]: 

= and jU { s m) =ju~ +RT In M 0 [4.10] 

M 0 is the molar mass of the pure solvent and x 0 its molar fraction, 
practically equal to 1 for dilute solutions. 

If, now, the unit of composition is the concentration (molarity in moles 
per liter), with v 0 denoting the volume of solvent, then if the solution is 
sufficiently dilute we find: 


c „ = - 


[4.11] 


and therefore the chemical potential of the solute becomes: 
M, = v: + R T In x s v 0 cy s U) = jU m> + R T In c/ s m) 


[4.12] 
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For the activity coefficient in reference (III), the molar solution: 

f s in) = f s n X^ riband nr=^\ [4-13] 

We shall also see in section 4.2.11 how to switch from using 
concentrations to using molalities. 


4.2. Debye and Hiickel’s electrostatic model 

The purpose of Debye and Hiickel’s electrostatic model is to give an 
account of the contribution of the effect of charges on the properties of 
solutions. When it was first published, this model was considered a model of 
an ionic solution in itself, but applicable only to sufficiently dilute solutions 
of strong electrolytes, which are fully dissociated - in practice, with ionic 
strengths less than 0.01 mole/1. A resurgence of interest in this model 
occurred when it was understood that the interactions due to the charges 
were not the only interactions at play in solutions, and that the model needed 
to be used in conjunction with another, integrating the molecular interactions 
over short distances. Debye and Hiickel’s model has thus become the 
standard way to express the contribution of the electrical charges to the 
properties of any solution containing ions. 

Numerous publications have given developments of Debye and Hiickel’s 
theory, some of them more condensed than others, and some more rigorous 
than others. In view of the importance of this model, which is used 
unanimously in all models of ionic solutions, we have chosen to base 
our presentation here on that given by Fowler and Guggenheim. Indeed, 
the method developed by these authors shows numerous advantages, 
including the fact that it is extremely rigorous in regard to the fundamental 
hypotheses - particularly hypothesis 4 (see section 4.2. 5. 4). This 
development also establishes criteria of self-consistency which must be 
respected by the accepted approximations. We have chosen to present this 
model by using the concentrations (molarities) in moles/liter to express the 
composition of the solutions. We shall see later on (in section 4.2.11) what 
becomes of the expressions derived from the model when the compositions 
are expressed in molalities (moles/kilo). In addition, if we use the SI units, 
the expressions of electrostatics will respect this system, in which the 
electrical permittivity of the medium is given by: 

£ = £ 0 D 


[4.14] 
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D is the relative permittivity or dielectric constant of the medium: it is a 
dimensionless property whose value is around 80, for example, in the case of 
water at ambient temperature. £ 0 is the electrical permittivity of a vacuum. 

Its value is 8.85 x 10' 11 SI, which equates to around 10“ 9 / 36;r . 


4.2.1. Presentation of the problem 

Thus, we consider an encompassing solution, in a continuum with 
dielectric constant D (which, as we shall see, will be that of the solvent). 
This solution contains 5 types of ions 1,2, ,.,s with positive or negative 
electro valences z l ,z 2 ,...,z i ,...,z s at the concentrations (molarities) 
clcl,...,cl...,c° s meaning that the volume V contains N l ,N 2 ,...,N i ,...,N s 
ions or n l ,n 2 ,...,n i ,...,n s ions-moles of each type. Each ion i carries the 
charge z,e, where e is the elementary charge, which is 1.602xl0" 19 
Coulombs. The electrical neutrality of the collection is expressed by any one 
of 

the following three equations: 

Y N i z i = 0 or Y n i z i = 0 or X c ° z , =0 [ 4 -15] 

s s i = 1 

The overall density of an ion j in the volume of the solution is 
N,V = N a c 9 / 1 000 (if the concentration is expressed in moles per liter). If 

we choose an elementary volume 8(0 at a point M of the solution in the 
vicinity of any ion k (Figure 4.1), then following thermal agitation, at all 
times, ions of the different species enter into that element, so that at a given 
time, there may be an excess of positive or negative charges. The local 
density of ions j around the point M will be SN . / 8 ( 0 = N a cj M) / 1 000 (where 

Cj is the concentration in molecules per liter around the point M). 

Nonetheless, on average, over time, the local concentration of an ion j is not 
equal to its overall concentration in the solution. Thus, if the point M is a 
neighbor of a cation, it is probable that the density of anions present will, on 
average, be greater at point M. Put differently, the electrostatic forces over a 
long distance lead to an ordered structure of the solution, the effect of which 
is that the concentration of a species at any given point is not equal to its 
overall concentration in the solution. We recognize the concept of local 
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composition, as described in Chapter 3. We write the following for the 
concentration of ions j at point M: 

c °j = c \ M) g kj ( r ki > z k z j T ) [ 4 - 1 6 ] 



The function g kj (r kJ ,z k z j T) is a distribution function which depends on 

the two species k and j, on the interactions between those species and 
particularly, on how far apart they are r kj , but also on the temperature, 

because thermal agitation works against the ordering of the solution. The 
distribution function will generally tend toward zero if the repulsion between 
the ions k and j is strong, whereas if there is a strong force of attraction 
between those two ions, this function will tend to have a high value. 

At a point M, the charge density would be: 

p = f j ^z i e = c[ M) Y J z >-c {M) Y J \ z -\e [4.17] 

*=1 i i 

We shall use the notation c[ M) and c[ M) to denote the local 
concentrations of positive and negative charges at that point M. 


4 . 2 . 2 . Notations 

To clarify our arguments, in this presentation of Debye and Huckel’s 
model, we shall use a particular system of notation, which we shall begin by 
presenting here. 
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Suppose we have a value Q that quantifies a property of the collection of 
ions and is dependent on their distribution within the volume, i.e. on the 
configuration of the system. Just as we have done up until now, we shall use 
(0 to represent the mean value of Q across all possible configurations. If w 
denotes the potential interaction energy of the system, then by virtue of the 
first mean value theorem, we shall have: 

(fi)f-Jexp (diy) A '=J...j0exp (d<y) A [4.18] 

V V J V \ / 

We use the notation (Q) a for the mean value of Q across all possible 

configurations around an ion a when that ion a is held still in a given 
position. Thus, we shall have: 

(2)J"\f ex P “wV ( d ®) A ~‘ = f-ffiexp (dfy) Y '‘ [4.19] 

v V ^-b-* J v V ^b-' y 

The integrations are then performed on the coordinates of all the ions 
other than the ion a. 

4 . 2 . 3 . Poisson ’s equation 

Hereinafter, we shall suppose that we can define, at each point of the 
solution, for any configuration at all, an electrostatic potential *F and a 
charge density p which obey Poisson’s equation. Thus, in light of relation 
[4.16], at a coordinate point r which is not occupied by an ion, we would 
have: 

V 2 W(r) = [4.20] 

e 0 D 

If we now average this expression over all configurations with the fixed 
ion k , in light of the notations introduced in section 4.2.2, we shall have: 

{PW) k 

£ 0 D 


V 2 <W(r)) t 


[4.21] 
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4.2.4. Electrical potential due to the ionic atmosphere 

We have agreed (expression [4.1]) that the interaction energy could be 
considered to be the sum of an energy due to the electrical charges 
and an energy due to the molecular interactions in proximity. We now 
suppose that, if the molecular interactions were the only forces at play, the 
solution would be accurately represented by the model of an ideal dilute 
solution (see section 3.1.3) and that the deviation from this model is due 
merely to the long-distance interactions w lCh> caused by the electrical 
charges. 

Because our solution is dilute, the energy 1111,1 1 is, in fact, due only to 
the molecular interactions (short distance) between molecules of solvent and 
ions of solute, and therefore is independent of the different relative 
configurations of the ions. 

The Helmholtz energy of the real solution with the charged ions would 
exceed that of the hypothetical solution, wherein the ions are not charged 
with an amount F (Ch) such that: 


r 

F <Ch) = -k B 7’lnj'...Jexp 

V 


w (Ch)N 

KT , 


n(d®,f'+iVk B rinr 


[4.22] 


in addition, the interaction energy due to the electrical charges can, itself, 
be broken down into two terms: 

- the energy vv lScl/ ’ due to the charging of each ion in the absence of all 
the other charges; 

- the energy w idec) due to the electrical interactions between the ions: 


w ( ch > = w < Self > + 


[4.23] 


The Helmholtz energy F iCh) can thus be written: 


p(Ch) _ (Self) + p(elec) 


[ 4 . 24 ] 
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We can therefore write the excess Helmholtz energy due only to the 
electrical interactions in the form: 


r 

p(elec) _ _k B rln|...| exp 

V 


w (dec) ' 


nid^f'+^VkBrinr 


[4.25] 


Let us derive relation [4.24] in relation to the charge of the ion k:e\z k |, 
after having used expression [4.25]. We have: 


dw (sm dF (e,ec) 

1 

ed z k ed z k 


r rdw ,Ch) 

J J ^ r ex P 


f W (ch) ^ 


ed 




(d®)* 


{...{exp 


f w <Ch)A 


V k B^ J 


(d®)* 


/ dw (Ch) 

\ed\z k 


[4.26] 


The value 



is the mean increase in energy of the ensemble of 


particles per unit charge of k. Thus, this mean is equal to the electrostatic 
potential at the point occupied by the ion k and therefore, in view of our 
system of notation (see section 4.2.2), we can replace relation [4.26] by: 


aw Self) 3F (elec) 

1 

ed z k ed z k 


(no)}„ 


[4.27] 


However, 3rv (Sclf) I ed\z k \ is the fraction of the electrostatic energy in k 
due to the ion k itself, and therefore this term is independent of the other 
ions. If we subtract that term from , we obtain the fraction of the 

mean electrostatic energy occupied by k and due to all the other ions. This is 
what is known as the potential due to the ionic atmosphere. It is written as 
^,,so: 


3w (Self) 

ed z k 


n = (^(0)). 


[4.28] 
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and therefore, we find the following for the contribution JF ' (dec) ; 

■Wetec) 

-^r-r = ^k [4-29] 

e °\ z k\ 

If we take into account the contribution of each ion i in the set of ions, we 
can always write: 


dF (elecJ =e£V' d|z,. 


[4.30] 


Each term F i is a function of all the charges z x , z 2 , . . . , z„ ... of all the 
ions. 

However, the Helmholtz energy is a function of state, so we can apply 
Schwartz’s equation, meaning that we have: 

a*/ 7 , a f , 

— ± = L T4 3 11 

a N d \ Z k 

This relation will be extremely useful because, as it is very general, it can 
be used to test the self-consistency of any approximate expression of the 
potentials F t . 


4 . 2 . 5 . Debye and HuckeVs hypotheses 

Debye and Hiickel’s model is an approximate evaluation of the potentials 
F ) and of p(r ) , in order to solve Poisson’s equation [4.24]. This model is 
founded on five hypotheses. 

4.2.5. 1. Hypothesis 1: shape of the ions and nature of the medium 

The ions are all considered to be hard spheres with the same radius a and 
different charges. The medium in which they are bathed is considered as a 
continuum with the dielectric constant D chosen as being that of the pure 
solvent, independently of the presence of the ions. 
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4.2. 5. 2. Hypothesis 2: pairwise interactions 

We suppose that the energy w (elec) , which is due to mutual electrostatic 
interactions between the ions, is the sum of all the potential energies of 
interaction of the pairs of ions which it is possible to envisage on the basis of 
the ions present in the solution. We write this condition as: 

w< e' eC, =IIX [4-32] 

k=l j*k 

If we take account of relation [4.14], the pairwise energies £ kj are such 
that: 


£ v = 


\ z A e 


4ne 0 Dr kj 


for r kj > a 


[4.33] 


and 


s kj =oo for r kJ <a 


[4.34] 


This last condition simply expressed the fact that the centers of two ions 
cannot come any closer to one another than a distance 2a. 

4.2. 5. 3. Hypothesis 3: Boltzmann distribution 

The third hypothesis is to accept that the distribution function of the ions, 
g kj {r kj ,z k z jT ) , in the space is such that at the distance r kj (greater than or 

equal to 2a) of an ion k, the concentration of ions j is given by: 


f 


c. = c , exp 


KTj 


AT f 

or SN, = — l Sojcx\) 

1 V 


KT; 


[4.35] 


This is the form of Bolt z mann’s distribution. 

NOTE 4.1.- The energy £ kj has multiple meanings. Indeed, it is: 

- the mean Helmholtz energy of an ion j at a distance r kj from an ion k; 
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- the work needed to bring an ion j from an infinite distance to the 
distance r kj from the ion k (this work is averaged over all the configurations 
of the other ions); 

- the potential energy, whose derivative gives the mean force exerted 
between two ions k and j, set a distance r kj apart. 

From this third meaning, we can deduce: 


£ kj - £ jk 


[4.36] 


This expression is not clear for an energy averaged over all the 
configurations. 

Based on the expressions [4.35] and [4.17], we can calculate the charge 
density at the distance r from the ion k. We can then easily calculate: 


(Plr))k = 2>. 


SN. 
; e - 

8m 


= wX-v v , cx p 


V 


j^k 


k B r 


[4.37] 


By substituting this back into Poisson’s relation [4.21], we obtain: 


vlW,) 


s 0 DV j* k 


Z z j N .\ ex p 


KT; 


[4.38] 


It is envisageable to solve this equation when we know the relation that 
exists between and £ kj . 

4.2. 5. 4. Hypothesis 4: relation between ^F(r)') and e kj 

Debye posits a fundamental approximation which is a relation between 
the mean energy £ kj and the mean potential l*F(r)V, and accepts that we 

have: 


e « = e \ z j\{' F lr))k 


[4.39] 
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Using equation [4.36], we can deduce that: 

=h\(nr))j [4.40] 

Irrespective of the pair of ions in question, this relation [4.41] can only be 
verified if: 


{nr)) k _(nr))j 


[4.41] 


If we feed relation [4.39] back into expression [4.38] we obtain Debye’s 
equation: 




e 0 DVf? k 


Zj N j exp 


g | z /|( y ( r )), 

k B J 


[4.42] 


To simplify the notations, we shall write F(r) instead of (F(r)) and 
rewrite relation [4.42] in the form: 


V 2 (^( r)) 


£ 0 DV 


ZhK ex p 


e z t 'Fir) ' 

KT y 


[4.43] 


The relation is often called the Poisson-Boltzmann relation. 

4.2. 5. 5. Hypothesis 5: primacy of thermal agitation 

Although equation [4.43] can be solved, Debye’s final hypothesis, the 
aim of which is to simplify this resolution, is to consider the case where the 
thermal agitation energy is much greater than the interaction energy: 

£ kJ «KT [4.44] 


On the basis of equation [4.39], the above equation can be expressed by 
the inequality: 


ezj F{r ) « k B T 


[4.45] 
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This approximation enables us to replace the exponential in relation 
[4.43] with a limited expansion, which Debye limits to the second term - 
i.e.: 


r 

exp 

V 


k J , 


e N( y(r) )* 

k B r 


ejzj^(r) 

kt 


[4.46] 


Indeed, we can easily see that the first term of the expansion disappears 
in expression [4.43], because of electrical neutrality, which imposes the 
following condition: 

2>^=0 [4.47] 

(=i 

Thus, we can write the Poisson-Boltzmann equation [4.43] can be written 
in the form: 

V 2 ('F(r)) = /c 2 'F(r) [4.48] 


k denotes a constant defined by: 


x- 2 = 


e 0 DVk B T 



N a e 2 

£ 0 DVk B T 



where k > 0 


[4.49] 


If we use the concentrations to measure the composition of the solution, 
and if these concentrations are expressed in mole/1, rather than in mole/m 3 , 
relation [4.49] is written as: 




n y 


1000 £ 0 Dk B T ~T } 


Z 0 2 

C,. z t 


[4.50] 


In this new definition of the constant 1 C , we see the appearance of the 
ionic strength / as defined by relation [A.2.49]. Thus, this constant can also 
be written as: 


k 2 =- 


2N,e 2 


1000f 0 Dk B r 


[4.51] 
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We shall now easily be able to integrate equation [4.48] by using the 
spherical coordinates. 


4.2.6. Debye and HuckeVs solution for the potential due to the ionic 
atmosphere 

In spherical coordinates and spherical symmetry, the Laplacian contains 
only the radial term, which is: 


V 2 (W(r)) 


1 



df(r)l 


r 2 d r 


[4.52] 


Thus, in light of equation [4.41], the equation to be solved is: 


d 2 W(r) | 2 d*F(r) 
dr 2 r dr 


= tc ll F(r) 


[4.53] 


Equation [4.53] is a second-order differential equation, which implies that 
the solution must satisfy boundary conditions to set the two integration 
constants. 


The first condition is that the potential Y / (r) and the field, i.e. the 


derivative of the potential 


dW(r) 

dr 


, become null as r tends toward infinity. 


The second condition is that the electrical induction be constant at the 
boundary between the inside and the outside of the ion, meaning that for 
r = a. 


The solution of equation [4.53], which remains finite as r tends toward 
infinity, will be of the form: 
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K is a constant determined by the continuity of induction. The radial 
component of the electrical induction for r>a, according to relation [4.54], 


9*F(r) 


Ane n DK 


(1- /rr)exp(-Kr) 


[4.55] 


However, in the vicinity of the boundary, within the ion, induction is due 
only to the charge on the ion, and therefore its value is: 


dT{ r) 


[4.56] 


Let us write the continuity of the induction for r = a , which is: 


Ane„DK , , , , 4 ne z k 

5 — (l-*»)exp(-x») = — — — 

a Ana 


[4.57] 


By drawing K from expression [4.57] and substituting it back into 
relation [4.54], we obtain the potential as a function of r. 


V (r) = exp [~K(r - a)] 

Ane 0 Dr(\ + m) 


[4.58] 


and in particular, for r = a, this solution is: 


na) = ^ r 

Ane 0 Da{\ + Ka) 


[4.59] 


If we now subtract from IF(a) the potential due to the ion k for r = a, we 
obtain the potential *F k due to the other ions, i.e. the potential due to the 
ionic atmosphere: 


e\z t \ k 


Ane 0 Dr(l + tea) Ant'JJa AneJ) 1 + m 


[4.60] 
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4.2.7. Charge and radius of the ionic atmosphere of an ion 

Consider a sphere centered on an ion k with radius r. The areic charge 
density at the surface of the sphere is calculated on the basis of the mean 
charge density: 

a s (r) = Anr 2 {a) k [4.61] 


By taking account of Poisson’s law and of expression [4.48], we find the 
following for the areic density: 


c 7 S (r) = 4 Kr 2 £ 0 DK 2 'F(r) 


[4.62] 


However, *F(r) is provided by expression [4.58], which gives us an areic 
charge density on the sphere of: 


a sW = ; 


rice I. 


1 + Ka 


- exp( Ka ) exp(— ter) 


[4.63] 


If we study this function of r, we see that it exhibits a maximum for: 
d[rexp(-xr)] 


dr 


- = 0 


[4.64] 


By solving equation [4.64], we find that this maximum is obtained by: 

[4.65] 


r A =l/K 


This radius r A is known as the radius of the ionic atmosphere around an 
ion k. It is the distance at which we find the maximum of the electrical 
density of the ionic atmosphere of the ion k. 

Figure 4.2 shows the shape of the areic distribution as a function of r for 
two values of the ionic strength (7=1 and I = 10' 2 moles/1). We can see that 
the maximum point of the curve shifts toward larger radii if the ionic 
strength decreases, and that the amplitude of the corresponding peak also 
decreases. 
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In water at 25°C, we calculate (in meters): 


30.8x10" 


r, =' 


V7 


[4.66] 


Table 4.1 gives the value of the radius of the ionic atmosphere for these 
two values of the ionic strength. 


Ionic strength (moles/1) 

Radius r/nm) 

0.01 

308 

1 

30.8 


Table 4.1. Radius of the ionic atmosphere for a few values of the ionic strength 


If, for the radius of the ion, we choose the value of 30 nm, this means that 
for an ionic strength of 1 mole/1, in water at ambient temperature, the radius 
of the ionic atmosphere is of the order of magnitude of the radius of the ion, 
so: 


for/= 1 mole/1 r A = — = a [4.67] 

K 

Note that, as k is proportional to the square root of V, the radius of the 
ionic atmosphere will increase with the square root of the volume and 
therefore of the dilution. 
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4 . 2 . 8 . Excess Helmholtz energy and excess Gibbs energy due to charges 

Considering condensed solutions, we ignore the variations in volume, and 
therefore we can treat the excess Gibbs energy due to the charges as the 
same thing as the corresponding Helmholtz energy, and returning to relation 
[4.30], we write: 


dG felec) =dF (e,ec) =e£¥' i d|z j 


[4.68] 


In order to integrate this equation, we shall imagine a virtual path which 
consists of progressively charging the ions. A point on that path would be 
characterized by a fractional extent of the charge A between 0 and 1 . At that 
point, the charge of the ion / would be Aejz.j, and we integrate relation 
[4.68] in relation to A : 

i 

^«(elec) fiw (A)e\z i \dA [4.69] 

0 i 


Of course, the envisaged path is purely a calculation ploy, because a 
charge less than the elementary charge is inconceivable. 

Integration is done in the following manner: 


G (eiec) = _ y N,zje-K r A-dA = _ Na y 
^ 4 Jie 0 D { 1 + A tea ^ 


ntfe 2 K 

- L - L t{ko) 

12 7T£ 0 D 


[4.70] 


The function r(x) is defined by: 




log(l + x)-x + y 


[4.71] 


which, whenx is smaller than 1, we can expand into: 
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Table 4.2 contains several values of the function t[x) . 


NOTE 4.2 - From the approximation which we have just added, we can see 
that the expression envisaged for the excess Gibbs energy is only acceptable 
for Ka « 1 . This means that, for a given radius (around 30 nm), our Gibbs 
energy is acceptable only for values of k such as 3x 10“ 10 ; (fc« 1) , which 
corresponds to an ionic strength less than 1 mole/1. However, this restriction 
is not a true restriction because, for such values of the ionic strength, it has 
certainly been a long time since Debye’s hypotheses have been proved to be 
unacceptable. 

Thus, we can keep only the first term in expansion [4.72] and write the 
following for the excess Gibbs energy: 


^r(elec) 


N.zfe K 


1 


=--y 

12 ke 0 D 3 i 


Njzy 

V £ o D 




[4.73] 


By switching to concentrations in moles per liter, we find: 


^r(elec) 


-N a v^ 

a /_( , - 


2 2 
n,z, e K 


\27C£ q D 




■ n , c^yv ^' 2 

' lOOO^D 


[4.74] 


By revealing the ionic strength /, we obtain the following for the 
electrostatic excess Gibbs energy: 


^r(elec) 


2 V 2 r 

v 


N a e 2 F7 

1000£- o D 


\ 3/2 

) 


[4.75] 


Thus, we shall use the formulae [4.75] and [4.76], drawn from the first 
form of relation [4.66]: 


^r(elec) 


1 N , d /cVe 2 
3 1000.4^f 0 Z) 


I 


0 2 
c, z. 


V > J 


2N . A V>ce 2 jj 
3000A7TE 0 D 


[4.76] 



138 Modeling of Liquid Phases 


4.2.9. Activity coefficients of the ions and mean activity coefficient of the 
solution 

By deriving equation [4.70] in relation to the amount of component k in 
the solution, we can calculate the excess partial molar Gibbs energy of the 
ion k due to electrical charges: 


In Ytu) = At = 


9G ( elec) 

dn,. 


2 2-v T 

N a 


Z N ^t 


2 2 

e 


87T£ 0 D 1 + Ka 12 7T£ 0 D 2V 


y m(k) 


a(tca) 


This gives us the activity coefficient of the ion k\ 


In 


.(//) _ . 




2 2 23 

Z7e K K 

a + 

%7t£ Q Dk B T 1 + Ka 24 7T£ 0 DV 


K,a^(Ka) 


(t(x) is the function defined by: 


C7(X) = - 


1 + jc — - 


1 + x 


- 2 ln(l + x) 


[4.77] 


[4.78] 


[4.79] 


or 


1-3 — x + 3— x 2 -3—x 3 +3— x 4 - .. [4.80] 

4 5 6 7 

A number of values of (j(x) are given in Table 4.2. We can see that at a 

very high level of dilution, if Ka « 1 the function a(Ka) tends toward 1 

and in the conditions we chose previously, the activity coefficient is such 
that: 


In i/ /7 > = - 

111 fk{c) 


2 2 

z k e 


K 


2 2 
z.e 


2N a eA 
lOOO£ 0 Dk B T 


&7t£k B T 1 + Ka 


2D£ 0 k B T 


1 + a , 


2 Ml 

mO£ 0 Dk B T 


[4.81] 
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This function is usually written in the form of the decimal logarithm: 


=~ Bz , 


41 

1 + Aa4l 


[4.82] 


X- 

X 

x(x) 

o(x) 

0.0000 

0.0000 

1.000 

1.000 

0.001 

0.0316 

0.976 

0.954 

0.002 

0.0447 

0.967 

0.936 

0.003 

0.0557 

0.960 

0.922 

0.004 

0.0633 

0.954 

0.912 

0.005 

0.0707 

0.949 

0.902 

0.006 

0.0775 

0.945 

0.893 

0.007 

0.0837 

0.941 

0.886 

0.008 

0.0894 

0.937 

0.879 

0.009 

0.0947 

0.934 

0.871 

0.010 

0.1000 

0.931 

0.866 

0.100 

0.3162 

0.811 

0.659 

0.200 

0.4472 

0.752 

0.569 


Table 4.2. Table of values of functions T(x) and ofx) 


with the following meanings, for the term B (sometimes called Debye and 
Huckel’ s constant) and the term A : 


]JlOOOe 0 Dk B r 

2De 0 k B T 


B = 2.303 


[4.83] 
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A = 


2 N„e 


1000£- OJ Dk B r 


[4.84] 


Let us recap the values of the electrostatic constants: 

6-0 = 36^ and e = 1 ’ 602 - 10 " 19C t 4 ' 85 ] 

The law [4.82] is called Debye and HiickeTs extended law (as opposed to 
the limit law, which we shall see in section 4.2. 13). 

At 25°C, by choosing to use the dielectric constant of water for aqueous 
solutions (D = 78.54), we calculate: 

5 = 0.511 1 1/2 mole -0 5 [4.86] 

A = 0.3287 xl0 10 l 05 mole _05 nf 1 if a is in meters. [4.87] 

NOTE 4.3.- It is worth noting that if we choose the value of the radius as 
a = 3.04xl0" 10 m , we find Aa ~ ll°' 5 mole~ 05 , and relation [4.82] is then 
rewritten as: 


lo g f«c) = ~ Bz , 


77 

i + 77 


[4.88] 


The “electrostatic” chemical potential of the solvent is calculated, for 
instance, as follows: 


(elec) _ 

Mo 


3G (elec) _ 3 G (elec) d/c dV 

3 n 0 3 k dV 3 n 0 

Using relation [4.77], we obtain: 

2>,a 2 < 


[4.89] 


2 2 

", e 


(elec) _ _j_ 

Mo , 


K 


\27C£ q D 2V 


v m oM Ka ) 


[4.90] 
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In the context of very dilute solutions, a(ica) ~ 1 , this potential becomes: 


./elec) _ 

M o 


lAfz:V 


K 


12 7T£ n D 2V 


[4.91] 


Vg is the molar volume of the pure solvent. The solvent’s activity 
coefficient can thus be calculated, and we shall verify that it satisfies the 
Gibbs-Duhem relation. 

We know that it is impossible to determine the coefficients of the 
individual ions by experimentation, but that we can determine the mean 
coefficient Y±(c) of the solution. In light of electric neutrality: 


v + z + +v_z_ = 0 [4.92] 

Using relation [A.2.45] for the mean activity coefficient, we can easily 
calculate: 


MS=- 


z , z e 


87T£ 0 Dk B T 1 + Ka 


[4.93] 


With D = 78. 54 for water at 25°C, we calculate the following, with a 
being in meters: 


log YJu) = — 0.511z + z_ 


V7 

1 + 0.3287 xlO -10 a V7 


[4.94] 


Although the dielectric constant varies with temperature, the product DT 
which is the only one that plays a role practically does not change between 
18 and 25°, and therefore the values of the coefficients A and B do not either. 


4 . 2 . 10 . Self-consistency of Debye and Hiickel’s model 

Any model of an ionic solution must satisfy two categories of self- 
consistency criteria: thermodynamic criteria and electrostatic criteria. 
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4.2. 10. 1 . Thermodynamic criteria 

We have already checked, at the end of section 4.2.9, that the Gibbs- 
Duhem relation was verified. With regard to the other criteria: 

- the model must satisfy the symmetry of the characteristic matrix: 


dn , (hr 

J 1 


[4.95] 


If we apply this relation to the expression of the chemical potential, given 
by relation [4.77], we note that this relation is respected, provided we take 
both terms from equation [4.77], because the first term on its own is not able 
to satisfy this criterion. Certainly, the second term may be numerically 
negligible, but it cannot be taken as equal to zero. Hence, Debye and 
Hiickel’s model satisfies this criterion. 

- the terms in the diagonal of the characteristic matrix must be positive: 

^>0 [4.96] 

3 / 7 ,. 


We can see, from relation [4.70], that Debye and Hiickel’s model satisfies 
this criterion. 

-the excess Gibbs energy must be homogeneous of degree 1 (like all 
extensive values) in relation to the amounts of material. 

If we examine relation [4.50], we can see that k is of degree Vi in 
relation to the concentrations and that \fl is also of degree Vi, and therefore, 
in relation [4.76], the “electrical” excess Gibbs energy is of degree 1 and the 
function is homogeneous because each quantity TV, is the product: 

TV,. = N a nx i [4.97] 

and therefore the quantity n can be factorized. 

4.2.10.2. Electrostatic criteria 

We have three electrostatic criteria to satisfy. 
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The first criterion needing to be respected is the proportionality between 
the potential due to the ionic atmosphere and the charge. Relation [4.58] 
shows that Debye and Hiickel’s solution satisfies this criterion. 


We have encountered the criterion given by relation [4.31] that must be 
satisfied by the potential of the ionic atmospheres of the ions. Let us now 
examine whether Debye and Hiickel’s solution conforms to this criterion. By 
writing: 


_9W a die 
dz j die dzj 


[4.98] 


Fowler and Guggenheim show that the derivative d*F k / dz j contains the 
product ZjZ k , and also only depends on z ; . and z k by way of the coefficient 
k , and therefore this derivative is symmetrical in relation to j and k , and 
thus the criterion [4.31] is verified. 

The third criterion is that the solution must respect electrical neutrality. 
To determine this, we calculate the overall charge of the atmosphere 
surrounding an ion k. We shall consider two spheres centered on the ion k, 
with respective radii r and r + dr. The electrical charge contained in the 
volume between these two spheres is: 

d^=cr s (r)dr [4.99] 

By integrating this equation between r = a and an infinite value of r, we 
obtain the charge of the whole of the ionic atmosphere surrounding the ion k, 
meaning that the electrical charge of the whole solution decreases by that of 
the ion k. Thus, we have: 


^=jV,(r)dr [4.100] 

a 


Returning to expression [4.63] for the areic charge, we calculate that this 
overall charge is: 
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Thus, the charge surrounding the ion k is precisely the opposite of the 
charge on that ion, which shows that Debye and Hiickel’s law respects the 
condition of electrical neutrality. 

We have gone into detail in our examination of the respect of these 
criteria, because we shall see that certain corrections, proposed to improve 
the model’s conformity to experience, yield solutions which no longer 
respect all these criteria. 


4.2.11. Switching from concentrations to molalities 

Up until now, we have chosen to quantify the composition of solutions in 
concentrations (molarities) expressed in moles/liter. Very often, users of 
Debye and Hiickel’s model use molality values M,. , which are expressed in 
moles per kilogram of solvent. Between these two values, if the solution is 
sufficiently dilute (which is the case in the domain of validity of Debye and 
Hiickel’s model) to enable us to treat the volume of the solution and that of 
the solvent as one and the same thing, we can write: 


m /A _ c 

1000 ' 


[4.102] 


The coefficient 1000 stems from the fact that the molality is given in 
moles/kg and the molarity in moles/1, rather than in moles/m 3 . 

If we work on the basis of the second expression of the “electrical” 
excess Gibbs energy given by relation [4.76], using expression [4.102], we 
obtain: 


^»(elec) 


1^\nV k B T 


/ 9 9 \ 3/2 

' r N a M ; z,VM 0 A 

v f 1000f 0 D 


[4.103] 


Thus, we define an ionic strength in terms of molality, which is 
expressed, similarly to definition [A.2.49], by: 
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Thus, the “electrical” excess Gibbs energy can be written: 


^»(elec) 


2V2 N a M 0 e 2 C 

3^4 xVk^T 1000 e 0 D 


[4.105] 


Similarly, on the basis of expression [4.50], using relations [4.102] and 
[4.104], we write the coefficient k 2 in the form: 


KT =- 


2N a/ y 2 

iooOf- 0 z)k B r 


[4.106] 


and for the activity coefficient li nk ed to the molality values, we find: 


log ^7= 

1 + Ap 0 a^jI M 


[4.107] 


For water, though, we apparently have p 0 = 1 kg/dm 3 at normal 
temperature, so relation [4.82] with the values of the coefficients B and A 
given by equations [4.86] and [4.87], is conserved, with the ionic strength 
relative to the concentrations being replaced by the ionic strength relative to 
the molality values Im- 

Thus, equation [4.82] applies indifferently, whether the compositions are 
expressed in molarity or molality. The ionic strength is expressed, in each 
case, with the same unit of composition. 

In addition, we can express the activity coefficient in any one of the 
conventions, but using relations [4.10] and [4.13], if the solution is 
sufficiently dilute, the values of the different coefficients become the same, 
because we find: 

yf n) » r ! n > « y) m) [4.108] 

Definitively, all the activity coefficients are given by relation [4.82] with, 
in an aqueous solution, the values of the coefficients B and A given by 
equations [4.86] and [4.87], and choosing / or I M depending on whether the 
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compositions are expressed in molarity (or concentration) or molality. The 
same is true for the mean activity coefficients, using relations [4.86] or 
[4.87]. Thus, we no longer make the distinction between the different 
activity coefficients, simply writing them as y i , and similarly for the mean 
activity coefficients, which we shall denote by y ± . 


4.2.12. Debye and Hiickel’s law: validity and comparison with 
experimental data 

Experience tells us that ionic solutions only behave like ideal dilute 
solutions ( y ± = 1) below a concentration of 10" 4 moles per liter, whereas for 
a molecular solution, such behavior is acceptable below 10 2 moles per liter. 

Figure 4.3 shows a comparison of the variations of the mean activity 
coefficient of magnesium chloride as a function of the ionic strength. The 
points are obtained experimentally and the downward curve is obtained 
by application of Debye and HiickeTs law. The figure demonstrates that 
the model begins to deviate from the experimental values long before 
the ionic strength reaches one. In particular, for most electrolytes, the 
real curve exhibits an extremum which Debye and HiickeTs model never 
shows. 



Figure 4.3. Variations of the mean activity coefficient of 
magnesium chloride with the ionic strength 
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In practice, we consider that Debye and Hiickel’s law, expressed by 
relation [4.82], and in particular using its form [4.88] in water, can only be 
used for solutions in which the ionic strength is no greater than 10‘ 2 mol/l, 
although in certain cases, ionic strengths of 1 O’ 1 mol/1 can also yield 
acceptable results. Certain authors prefer to adjust the law [4.82] by 
adjusting two experimental conditions: 

- the extrapolation for an ionic strength of zero for which the mean 
activity coefficient is one; 

- the radius of the ion a with the experimental curve. 

By this method, we are able to apply the formula [4.82] to a broader 
range of ionic strengths. However, the values thus determined for a still need 
to be acceptable in terms of the physical meaning. Indeed, certain 
experiments yield very low values of the radius a - sometimes zero, and 
sometimes even negative values are required. The true value of a must be no 
lower than the ionic radius determined in a crystal of a corresponding salt. 

Nevertheless, let us remember that Debye and Huckel’s law supposes the 
ions are spherical, which can be accepted for sufficiently dilute solutions in 
which the radius of the ionic atmosphere is larger than a and therefore in 
which the ions are relatively far removed from one another. If the ions come 
closer together, the hypothesis of sphericity becomes trickier to accept for a 
large number of types of ions - particularly for polyatomic ions. 


4.2.13. Debye and Huckel’s limit law 

For very significant degrees of dilution (10" 4 </< 10' 3 ), the denominator 
term in law [4.82] becomes much smaller than 1, and we obtain what is 
known as Debye and Hiickel’s limit law. 

This limit law is written as: 


log y. =-Bzfyfl [4.109] 

The curve giving the mean activity coefficient, calculated using this limit 
law, is tangential to the origin of the complete curve given by equation 
[4.82], as illustrated by Figure 4.4, and can be used for the extrapolation of 
the real curve at infinite dilution. 
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Figure 4.4. Comparison of the complete law and Debye and Hiickel’s limit law 

The fact that this limit law can be obtained numerically by making a = 0 
in expression [4.82] does not mean that it must be attributed to an 
approximation in which the ions are assimilated to single-point charges. This 
approximation, which is sometimes presented, is totally erroneous, because 
it is easy to show that a collection of electrically-charged points is 
completely unstable; indeed, it is the existence of a non-null radius a which 
lends stability, because two ions can never come together completely, as is 
demonstrated by the second form [4.34] of the potential function. 

4.2.14. Extensions of Debye and Hiickel’s law 

The limits exhibited by Debye and Hiickel’s law in its application to real- 
world cases have led many researchers to attempt to improve it. These 
improvements were first made in the strict framework of the law [4.82], 
before branching out in a new direction, integrating the molecular 
interactions that we shall see in sections 4.3 and 4.4. 

One of the earliest methods involved adding one or more adjustable 
parameters to the law [4.82]. Thus, Hiickel, taking account of a variation in 
the medium’s dielectric constant with changing ionic strength, was led to 
propose the following relation: 



[4.110] 
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This relation is able to represent the minimum of the curve in Figure 4.3. 
Unfortunately, the term K cannot be calculated a priori, and must be 
determined experimentally, which takes away from the advantage of this 
type of formula. 


Along the same lines, after numerous measurements, Davies proposed to 
determine the mean activity coefficient of aqueous solutions (at standard 
temperature) using the following relation: 


log y ± = -0.5 1 1 z + z. 



0.3/ 


[4.111] 


Fie proposed a second version of the relation: 


log 7+ = ~ B ,n Z + z - 


1 + C 1 A m a 'Jh, 


+ C 2 


[4.112] 


Fie drew up lists of values of the constants C \ and C 2 for each family of 
electrolytes. 

Figure 4.3 gives the result obtained by that model for magnesium 
chloride. We can see that whilst the model does show the minimum of the 
curve, it quickly deviates from the experimental results and, therefore, is 
hardly any better than Debye and FluckeFs original model. 

A second way of improving Debye and FluckeFs model was to review 
some of its hypotheses, making them less stringent. 

An initial attempt was made by modifying hypothesis 5 (see 
section 4.2. 5. 5), increasing the number of terms chosen in the serial 
expansion of expression [4.46] or indeed by solving the complete Poisson- 
Boltzmann equation with the exponential term, which is possible. Of course, 
the expressions obtained are complex, but they all manifest a major 
shortcoming: they no longer satisfy the electrostatic criterion represented by 
relation [4.41], which greatly takes away from their interest. 

A second way to improve was to keep certain terms (at least two) in the 
expansion [4.72], but we have already pointed out in section 4.2.8 that the 
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degree of accuracy obtained became illusory in relation to the other 
hypotheses. 

A third way of improving the law was attempted by Bjerrum, who 
introduced interactions between the ions in the form of associated species 
between ions that are sufficiently close. We can also introduce terms linked 
to the shapes of the ions by introducing parameters resulting from the 
associations between the anions and cations present in the solution. Using 
this technique, Guggenheim was able to develop a law with multiple 
parameters, which can be deduced from the experimental behaviors of pure 
salts. 

It is clear that, ultimately, these latter attempts consist of taking account 
of the interactions between the ions over short distances. This observation 
has led to a new general way of modeling ionic solutions: now we no longer 
modify the law [4.82] or its derivatives, but instead supeipose on the effects 
of the electrostatic interactions, represented by that law, the effects due to 
the molecular interactions over short distances, as found in non-ionic 
solutions. It is this method which we are going to discuss in sections 4.3 and 
4.4 below - a method which enables us to move away from just strong 
electrolytes, to take account of mixtures of ions and neutral molecules, and 
therefore weak electrolytes. 


4.3. Pitzer’s model 

Pitzer’s model is a semi-empirical model based on Debye and Huckel’s 
extended law [4.82], where additional terms have been introduced in order to 
take account of the effects of the ionic strength on the binary interactions 
over a short distance. 

Although we noted that this operation rendered the model non-self- 
consistent, in order to improve its conformity to experimental results, Pitzer 
proposed an expansion of the exponential of relation [4.43], limited to the 
first three terms instead of two, as Debye and Hiickel did in relation [4.46], 
so: 


r 

exp 

V 


w | 


k J) 



k B r 




v.k B ^y 


[4.113] 
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If Ay and fl ijk respectively denote the terms describing the pairwise 

interactions between the species i and j and the three-way interactions 
between the species i, j and k, expressed in kg 2 .mole' 2 and kg 3 mole' 3 , the 
excess Gibbs energy is written as: 


G xs 

R T 


^r(elec) 


R T 


+ZZ W 1 , +E£IX* M < M A 


[4.114] 


^*(elec) 

is the term expressing the long-distance electrostatic interactions 

RT 

between the ions. It is a function of the temperature and of the ionic strength. 
This function stems directly from Debye and Hiickel’s hard-sphere model 
(see section 4.2). With the new hypotheses chosen by Pitzer, it is expressed 
thus as a function of the ionic strength: 


^f(elec) 

R T 


= -x 0 M 0 B m ^-]n(l + by[l) 


[4.115] 


B m is Debye and Hiickel’s constant, expressed in relation to the molality 
values and given, as is shown by the comparison of relations [4.82] and 
[4.107], by: 

B m =Byfc [4.116] 

The parameter B is given by relation [4.86]. The term b is an adjustable 
parameter which has been optimized and taken as equal to 1.5kg 1 2 . mole 12 
for all temperatures and all solutes. It is theoretically linked to the distance 
beyond which the forces of repulsion between the ions become significant. 
p 0 is expressed in kg/dm 3 . 

NOTE 4.4.- In ionic solutions containing molecular species, the parameters 
Ay must include the molecule-molecule interactions of the ions, excluding 

the effect of charge, and molecule-ion interactions. The parameter jU jjk must 

include the molecule-molecule-molecule interactions (again including the 
ions but excluding the effect of charge). 
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We consider that the solution contains n c types of cations numbered 
1, 2, ...c,... n c and n a types of anions numbered 1, 2, ...a, ... n a . We shall 
rewrite relation [4.1 14] to separately introduce the terms li nk ed to cations 
and anions. The excess Gibbs energy becomes: 



4x o M o S '" / ln( 1 + ZiV7) + 2 IZ M M 

Bca 

+ (lz c M 

+Z2X M c 

c c' 

( Pcc: +] -Y. H aVcc-a 

^ a 


a a' 

<Paa’+ ^Y^cVcau' 

^ c 



[4.117] 


This relation exhibits the advantage of introducing adjustable parameters 
with experimental data: B ca , C ca , (p cc . and (p aa .. 


The terms B ca are given by: 

B ca =^ 0) + ^ , /(^: ) V7) + A 2, /(^ 2) V7) [4.118] 

The coefficients /L' 1 ' , p^J and p ]] ' are assimilated to second-order virial 
coefficients, expressing the short-distance effects between the ions c and a. 
They are expressed in kg/mole. The function /(v) is defined by: 

2fl-(l + y)exp(-y)l 

/b) = -^ 2 [4-119] 

y 

The second-order virial coefficients have the following properties: 

[4-120] 

| Py = P n = 0 if i and j have the same sign [4.121] 


NOTE 4.5.- The London forces of repulsion prevent interactions between 
ions over a short distance. 
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The parameters C ca are given by the expression: 


C 


C 


[4.122] 


The parameters (p cc , and (p aa , are also functions of the ionic strengths 
given by: 

Vcc' = Vc’c=0 ( £ + fi2V) [4.123] 

and 

Vaa'=Va'a=ffi+ffiV) M 24 ] 

The coefficients fifj and j8 ( c l J are functions of the ionic strength and of 
the charges on the ions. We have seen that these coefficients are null if they 
refer to cations or anions of the same charge (z a = z a or z c = z c ). 

Thus, the new form [4.117] of the excess Gibbs energy brings into play 
various parameters that are adjustable to each ionic solution studied. These 
parameters may be classified into two categories: 

- so-called “binary” parameters. This category includes the parameters: 

and and a , C* ca , which characterize binary solutions 

(with an anion a and a cation c in addition to the solvent); 

- so-called “ternary” parameters, which include the parameters and 
B ( V , w ,and v/ , . 

r^cc 9 t caa r cc a 

NOTE 4.6.- the parameters fifj and fi\'\ , which relate to pairwise 
interactions, are described as “ternary” because they characterize ternary 
solutions with an anion or a cation in common, as the electrolytes: mixtures 
of ac and a ’c or of ac and ac 

Pitzer’s model, as we have just described it, can be used for solutions 
whose ionic strength can be up to 6M. 
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For lower ionic strengths, it is possible to use simplified forms. In 
particular, we can often simplify the model by setting = 0 for any 
couple ca. 

The coefficients c/'J are such that, frequently, we can use the following 
values: 

- a ( c 'J -2 for all the electrolytes ca except those of the type 2-2(z c and 
z a , which are different from 2); 

- (/.[I ' =12 for all the electrolytes ca of the type 2-2 ( |z a | = z c =2). 

Therefore, in the highly-simplified version, there are now only three 
types of parameters to identify - namely /T' 1 ' , /?']' and C* a . In these 
circumstances, relation [4. 118] can be simplified to: 

- if al'J -2 - i.e. for all electrolytes ca except those of type 2-2: 


B = B' 0) -B m 

ca I ca r' ca 


l-3exp(2\/7j 


21 


[4.125] 


- if a t j - 12 - i.e. for all electrolytes ca of type 2-2: 


B =B' 0) -B" 

ca / ca # ca 


1 — 13exp(— 12>/7) 


111 


[4.126] 


For example, for a solution with only one electrolyte (a cationic species 
and an anionic species), it is easy to see that the excess Gibbs energy 
assumes the simplified form: 


4BJ 


R T 


hr(l + W7) + 2M a M c [S ca +z c M c Cj 


[4.127] 


The parameter C ca conforms to definition [4.125], and the parameter 
B ca , depending on the type of electrolyte, conforms to one or other of 
relations [4.125] or [4.126]. 
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The activity coefficients of the ions can be calculated in the following 
way: 


RTln^ 7 ' = — 

r ' dn t 


dG xs 9M, 

3M ; 


1 dG xs 
M 0 n 0 9M, 


[4.128] 


no is the amount of solvent of molar mass M 0 . 

Additionally, Pitzer’s model has been extended by Beuter and Renon to 
apply to solutions containing dissolved molecular species stemming from 
gases. 


4,4. UNIQUAC model extended to ionic solutions 

The extension of the UNIQUAC model to ionic solutions was first 
envisaged by Sander et al. in 1986. The basic idea is to think of the excess 
Gibbs energy of an ionic solution as being the sum of two contributions: 

-the term (j- a<e,ec > drawn from Debye and HiickeTs model, to take 
account of the long-distance electrostatic interactions between the ions (see 
section 4.2); 

- the UNIQUAC term as it is used for molecular solutions (see 
section 3.5), in order to take account of the short-distance interactions 
between molecules, between molecules and ions and between ions. 

Thus, from the excess Gibbs energy of the solution, we write: 

Qxs(ion) _j_ qxs(UNI) j-^ ^ 29 ] 

For a given species i, this relation results in a relation between the 
activity coefficients, given by: 

=\n^ ion) +\nrf II * mi) [4.130] 

Normally, the UNIQUAC model leads to the activity coefficients with the 
pure-substance reference ln^ (WV/) . To switch to the activity coefficient for 
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reference II - the infmitely-dilute solution - we use relation [4.2], so 
equation [4.130] is reformulated as: 

ln ^") = l„^)(.w + ln ^__ [4,131] 

The term ^ n '" NI) is the activity coefficient of the species i in infinite 
dilution. Thus, it is calculated by the UNIQUAC model by taking a molar 
fraction of one unit for the solvent (and therefore a molar fraction of zero for 
all the other components), so that: 

In = pm In y { j ,)iUN ' ) [4.132] 

* 0 ->l 


Thus, definitively, the parameters of the model are: 

- the parameters B and A (often taken as equal to 1.5) used by Debye and 
Hiickel as a function of the temperature; 

- the parameters linked to the UNIQUAC model, i.e. the structural 
parameters r, and q, and the energy parameters a u and a,, for each pair of 
components. 

The number of parameters needed for this model, therefore, is less than 
that needed for Pitzer’s model. 

The problem with electrolytes is that of the nature of the species that are 
actually present in the solution, which may be very different to those 
introduced when making up the solution. Indeed, we may see more or less 
complete dissociations, reactions with the solvent (water, for example), 
solvation of ions as introduced by Lu and Maurer (1993). These equilibrium 
states introduce new relations between the activity coefficients. For example, 
for the solvation of a cation c, the equilibrium is written as: 

/7cH 2 0 + c = c(H 2 0) ;;f 

This equilibrium introduces the expression of the law of mass action, 
which is the relation: 

_ A hc / he 


K 


[ 4 . 133 ] 
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In such cases, only an iterative method can be used to calculate the 
activity coefficients because, in order to find the molar fractions, we must 
already know the activity coefficients to use relations such as [4.133]. 

The UN1QUAC part of the model can be extended to the UNIFAC 
version in order to enhance the model’s predictive capability. 

This type of model can obviously be used for ionic solutions “with no 
solvent” such as pure salts in the molten state. 




5 


Determination of the Activity of a 
Component of a Solution 


In this chapter, we shall examine the methods for determining the 
activities, or the activity coefficients of the components of a condensed 
liquid solution whose composition is known. Note in passing that this 
determination of the activity may soon yield that of the chemical potential of 
the same component, if we know the chemical potential of that component in 
the reference solution. 

Remember that the activity and the activity coefficient of component of a 
solution depend simultaneously on the composition of the solution and its 
temperature. Consequently, a determination of the activity of a component is 
valid only at a given temperature and for a known composition of the 
solution. 

Remember, also, that there are several activities for the same component 
depending on the convention chosen to define it (see section A.2.5), namely: 

- convention (I) - pure-substance reference; 

- convention (II) - infmitely-dilute solution reference; 

- convention (III) - molar solution reference for all solutes. 

Hence, if it is not universal, each method must specify the convention in 
the context of which the determination is performed. 


Modeling of Liquid Phases, First Edition. Michel Soustelle. 

© ISTE Ltd 2015. Published by ISTE Ltd and John Wiley & Sons, Inc. 
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The methods for determining the activities (or the activity coefficients) 
can be divided into two main categories: 

- experimental methods, which are our focus in this chapter; 

- methods using a model of a solution - methods which stem directly 
from the study of different models (see Chapters 2, 3 and 4). 

The experimental methods which we are about to examine are based on 
the properties of the solutions - specifically on their behavior in a system at 
physicochemical equilibrium. The component whose activity we wish to 
measure is involved in a physical, chemical or electrochemical equilibrium, 
and the only u nk nown activity is that for which we are searching. Certain 
methods apply more to molecular solutions than to ionic solutions, whereas 
others are reserved for conductive solutions (ionic solutions or electrical 
conductors), and others are valid for all types of solutions. 


5.1. Calculation of an activity coefficient when we know other 
coefficients 

Before discussing experimental methods per se for determining the 
activity, we shall examine two methods which are, in fact, calculation 
methods based on knowledge of either the activity of the other species in the 
solution or the activity of the species in question at a different temperature to 
the desired conditions. 


5.1.1. Calculation of the activity of a component when we know that of the 
other components in the solution 

Suppose that, for all compositions between a known state (generally the 
reference state) and the composition under examination, we know the 
activities (or activity coefficients) of all the components of a solution in the 
same frame of reference, except for one of them, and we are going to 
calculate the u nk nown activity of that component for the chosen 
composition. This method is founded on the Gibbs-Duhem relation, which is 
valid regardless of the convention adopted. Thus, we obtain the activity (or 
the activity coefficient) in the chosen convention for the known values. 
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For the activity coefficients, at constant temperature and pressure and for 
a solution with N components (see equation A.2.23), this relation is written 
as: 


i> ; .dln^.=0 [5-1] 

;=i 

In light of relation [5.1], if we wish to determine the activity coefficient y 
of the component j whose molar fraction is Xj, we shall write: 

din y j =-Y J —^y k [5-2] 

k*J X J 

This relation is integrated between two states p and q. Each state is 
characterized by a composition of the solution: 




[5.3] 


In general, the lower bound p is chosen as the reference state, which gives 
a value of 1 for the activity coefficient /. , and the upper bound q is taken as 

the solution in which we are seeking the coefficient. In convention (I), the 
pure-substance reference, this gives us: 

din fP = -Y P— dln^ 7 ' [5.4] 

J 1 J ' ■ Jl v 

k*j A j 

Integration yields: 

In y~P = -Y f ? — din y i k I) [5.5] 

t'r Xj 


In convention (11), for the solvent, in the infmitely-dilute solution 
reference (which is to say, pure solvent), the relation is identical to the 
previous one: 


in>r=-zf;^ iw: n 

k* 0 ^0 


[5.6] 
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In convention (II), for the solvent, in the infmitely-dilute solution 
reference, the relation becomes: 

[5.7] 

k*s 

Expressions [5.5], [5.6] and [5.7] are usually used by numerical 
integration. 


5.1.2. Determination of the activity of a component at one temperature if 
we know its activity at another temperature 

If we know an activity coefficient at a certain temperature, it is possible 
to obtain the value for another temperature, knowing the dissolution enthalpy 
of the component as a function of the composition and the temperature. 

Let us employ the terms H* and H* respectively to denote the partial 
molar enthalpies of the component i in the reference solution and in the 
solution under study, for which the molar fraction of the component i is x. 
The application of Hel mh oltz’s second relation to the chemical potential of a 
component i (which is the partial molar Gibbs energy of that component) is 
written: 


d( Ml /T) _ H' - H* _ p d In yV 
dT T 2 dT 


[5.8] 


From this, we deduce the variation of the activity coefficient of that 
component: 


din yp _ /// - //; 
dT - RT 2 


[5.9] 


If we consider a solvent, regardless of the reference state, or a solute in 
the pure-substance reference (1), the above relation is written: 

d\nf n 
dT ~ RT 2 


[5.10] 
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The difference It - // appears as the opposite of the partial molar 
enthalpy of dissolution A diss H* of substance i from the pure state to the state 
of composition x. 

If, for a solute, we choose reference (II) - the infmitely-dilute solution - 
then relation [5.9] is written as: 


din _ Hr-H? 
dT R T 2 


[5.11] 


The numerator in the right-hand fraction can also be written as: 

Hr-H?=(W-h?)-(H?-h?) = A diss H~ - A diss H* [5.12] 


Thus, we find the difference between the variation of the partial molar 
enthalpy of dissolution of the pure substance in an infmitely-dilute solution 
and that of dissolution of the same pure substance in a solution of the chosen 
composition. These two values are given (Figure 5.1) on the curve 
representing the enthalpy of dissolution of the component as a function of its 
molar fraction. The first value A , liss H°° is the slope of the tangent to the 

origin of the curve, and the second value A diss H* is the slope of the tangent 
to the same curve at abscissa point x. 

It is possible to integrate the different expressions [5.1 1] and [5.12] if we 
know the variations of the enthalpies of dissolution with temperature. 

Usually, these enthalpies are considered to be constants within a 
reasonable range of temperature, so between two temperatures T and T, we 
can integrate relation [5.9] in the different cases in the form: 


7i(T) l n Yi(T') 


H* - H* ( 1 P 
R [t' T , 


[5.13] 


The partial molar enthalpies H* and H* or the enthalpies of 
dissolution A diss H* and A diss H°° are therefore constants, and one of the 
values of the activity coefficient jj** is known in the same composition. 
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Figure 5.1. Enthalpy of dissolution of a species 
as a function of the molar fraction 


5.2. Determination of the activity on the basis of the measured vapor 
pressure 

The vapor pressure methods derive directly from the properties of the 
equilibrium between the condensed solution and the vapor - a consequence 
of the equality of the chemical potentials of the component in question in the 
solution and in the gaseous phase at equilibrium with it. This equality gives 
us a very general relation between fugacity coefficient of component i in the 
gaseous phase @\ G) , the total pressure, the activity a\ L) of the component in 
the solution and an equilibrium constant K ( : LV) : 


c/y (, > p 

„(£) 


= K 


(LV) 


[5.14] 


In experimental terms, these methods only require measurements of the 
pressure and the composition of the gaseous phase. 

We can distinguish two approaches: 

- one, known as the direct method, which is primarily used with a pure- 
substance reference convention, meaning either for a component in a 
solution, for which we choose convention (I), or for the solvent of a solution, 
for which we choose convention (11); 
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-the other, which is based on the measurement of Henry’s constant, and 
is mainly used for solutes with the choice of convention (11). 


5.2.1. Measurement by the direct method 

For this method, we choose to adopt convention (1). Thus, the equilibrium 
constant for the liquid/vapor equilibrium is the same as the saturating vapor 
pressure of the pure liquid in question Ff , which is expressed by: 


0 iG) p 

Z_i £_ _ k (lg)(,d _ p o 


[5.15] 


In order for us to use this method, we need the pure component i and the 
solution to be in the same physical state so as to preserve the same standard 
state. Such is often the case, for example, with aqueous solutions of salts in 
which the water and the solution are in the same state, which cannot be said 
for the salt. Particular caution must be exercised with solid solutions: the 
solution and the solid i under examination must crystallize in the same 
crystalline system. 

The simplest case of the use of the direct method is found when the 
component being studied is much more volatile than the others, so that it 
alone practically constitutes the gaseous phase, which is then said to be pure. 
The conditions then allow us to treat fugacity and pressure as one and the 
same thing. Relation [5.15] then directly gives us: 


yLW) __ 


,(L)pO 


[5.16] 


The pressure P is often not very different from the saturating vapor 
pressure of the pure substance Ff . To improve the accuracy, therefore, it is 
preferable to work with a differential measurement. We use a pressure- 
differential sensor, which measures the difference in pressure that 
exists between two tanks placed in a thermostatic chamber at the chosen 
temperature. In one tank, equilibrium is established between the pure 
component and the vapor phase, and in the other, that equilibrium is 
established between the solution and the vapor phase. The saturating vapor 
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pressure of the pure gas is measured using an absolute pressure sensor in 
relation to the first chamber, or read from thermodynamic tables. 


When two components are notably present in the gaseous phase but it can 
be considered to exhibit perfect behavior both as a gas and as a solution, 
relation [5.15] gives us: 


y(L)U) _ 




JL) pO 


[5.17] 


In addition to the above measurements, the method requires us to know 
the composition of the gaseous phase, which can be found using a mass 
spectrometer. 

More generally, we can use the Lewis relation to characterize the fugacity 
of the component i in the gaseous phase on the basis of its fugacity 
coefficient in the pure state &‘ n G ' . That relation is: 

0(G) = 0°tO) x (G) 


Relation [5.15] then becomes: 


yU)(I) _ 


00(G) (G)p 


C W 


[5.18] 


The fugacity coefficient of the pure gas is supposed to be known or 
determined separately. 


5 . 2 . 2 . Method using the vaporization constant in reference II 

If we adopt convention (II) for a solute, the expression of the liquid/vapor 
equilibrium for the component i is written as follows at a given temperature: 


0 lG] p 

~/LXII) (L) 
r i x i 


= K UV) 


[5.19] 


If the gaseous mixture can be considered to be perfect, the equilibrium 
constant is the product of the saturating vapor pressure by Henry’s constant 
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(K^) of the component i, and relation [5.19] gives us the following for 
equilibrium: 


Y, 


= k (LV) =p°k 


[5.20] 


If the solution were perfect, it would obey Henry’s law, and would 
therefore satisfy the equation: 


Jr 


[5.21] 


This means that the constant K m is the limit, as the context x\ L) 
toward zero - in other words, if is equal to 1, with the ratio 


pV z) 


tends 
, so: 


K 


iH 


lim p 


[5.22] 


We determine the value of the ratio 0 ' at various decreasing values 

P x i 

of x\ L) , and extrapolate the curve obtained at x''/' 1 =0. The ordinate value 
than gives us the constant K m . 


The extrapolation is illustrated in Figure 5.2. The slope of the curve must 

P 

be horizontal. Otherwise, the measurements of the ratio - ' , have not 

r i X i 

been taken for sufficiently low values of the molar fraction of i in the liquid 
and thus found the validity of the limit law with a sufficient degree of 
accuracy. 


We can then feed that extrapolated value back into equation [5.20] and 
deduce the activity coefficient of the component i in convention (II): 


yUX H) _ 



,(G)i 


x\ L) P:K :f 


[5.23] 
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PJP, 


0 V (£) 



Figure 5.2. Determination of Henry’s constant 


As before, this law is easier to use when the component i is practically 
pure in the gaseous phase (x' G) =l), which is rare. If not, then a 
measurement of the composition of the gaseous phase is needed. 


5.3. Measurement of the activity of the solvent of the basis of the 
colligative properties 

A colligative property is a property of a solvent which depends solely on 
its molecular constitution rather than on the nature of the solutes in that 
solution. The existence of colligative properties is the consequence of the 
major dilution of the solute by the solvent. 

Thus, these methods are used to determine the activities of solvents. 
Hence, we always use the pure-solvent reference, irrespective of the 
convention chosen for the solution. Therefore, this convention no longer 
needs to be explicitly stated. 

Note that the method using the vapor pressures (see section 5.2.1) is 
already a colligative method. 


5.3.1. Use of measuring of the depression of the boiling point - 
ebullioscopy 

We suppose that the solute has a very low vapor pressure in comparison 
to that of the solvent. We can consider that the gaseous phase is pure and 
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contains only the vapor of the solvent. In the presence of the solute, the 
solvent’s boiling point is lowered. To calculate this depression, we write the 
equilibrium, which is expressed by the equality of the chemical potentials of 
the solvent in both phases, with the gaseous phase being under partial 
pressure of 1 atmosphere from the solution (at boiling point). Thus, we find 
the expression: 

d\ny 0 x ( 0 L) _ A„/? 0 ° 

dT R T 2 

To begin with, we can suppose that the vaporization enthalpy A v h° 0 is 
independent of the temperature. This being the case, it is easy to integrate 
this expression between the boiling point of the pure solvent T 0{Boil) and the 

boiling point of the solution T. If we use A7[ /Joi/) to denote the difference 
between T - T 0{Boil) , we find: 


In To = 


A v Ao A7 \ Boi i) 

R T 2 

S ^ 1 O(Boil) 


In 


1 - 1 - 


XL) 


[5.25] 


Thus, this relation enables us to determine the activity coefficient of the 
solvent at temperature T, which is not hugely different from T f)(fjoil) . 


If the solution is sufficiently dilute, the sum of the molar fractions of the 
solutes is much less than 1 and we can content ourselves with the first term 
of the expansion of the logarithm, and write: 


In To = 


A v h 0 AT {B oii) 

R T 2 

-^0 (Boil) 



[5.26] 


Now, if the enthalpy of vaporization depends on the temperature, it is 
sufficient to take it into account by integrating relation [5.24]. Note that this 
precaution is actually very rarely necessary, because the boiling points of the 
solvent and the solution are only a few degrees apart. 

For the moment, we shall restrict ourselves to applying relation [5.26] to 
molecular solutions. We shall discuss the application to ionic solutions in the 
next section (see section 5.3.2), with regard to the cryoscopic method. 



170 Modeling of Liquid Phases 


5 . 3 . 2 . Use of measuring of the depression of the freezing point - cryoscopy 

In a very similar manner to ebullioscopy, which we have just discussed, 
we can determine the activity, or the activity coefficient, of a solvent by 
cryoscopy, i.e. by looking at the depression of the freezing point of the 
solvent A T (P) owing to the presence of the solutes. If the solvent is pure in 

the solid phase, in the same conditions and with the same hypotheses, we 
establish a relation similar to relation [5.25] which, in the case of cryoscopy, 
is written as: 


In To = 


A fK AT (F) 

R T 2 

S ^ 1 0(F) 


In 


1-1-' 


XL) 


[5.27] 


Thus, we obtain the activity coefficient at a temperature T which is not 
hugely different from the melting point of the pure solvent T 0(F) . 


As before, if the solution is sufficiently dilute, the sum of the molar 
fractions of the solutes is smaller than 1, and we can content ourselves with 
the first term in the expansion of the logarithm and write: 


In To = 


A F h 0 A f F ' ) 

R T 2 

I ^ 1 0(F) 



[5.28] 


Thi s relation is easily applicable for molecular solutions. 

In the case of ionic solutions, a factor called the Van ’t Hoff factor , due to 
the dissociation of the electrolyte, comes into play in the sum of the molar 
fractions of the solutes. Thus, the situation is complex for weak electrolytes 
with incomplete dissociation, for which the dissociation coefficient is 
unknown. In order to resolve this impasse, it is possible to perform 
measurements using both ebullioscopy and cryometry at once. 

Let us take the example of a molecule of solute 5 which dissociates, 
giving rise to v s ions, and let a s be its degree of ionization. If we initially 
consider n s moles of non-dissociated solute, the total amount of solute n t(s ) 
after partial dissociation will be the sum of the amount of non-dissociated 
species and that of the ions: 


\ s) = n,(l-a,) + a s n s v s = n s [l + (v , -1 )a s ] 


[5.29] 
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In this case, the sum of the molar fractions of the solutes will be: 


y x \ ".[i + O-WR] 

. ,«o + « s [l + K -!)«,] 


[5.30] 


If the solution is sufficiently dilute, the second term in the denominator of 
the above fraction is smaller than n 0 (the amount of solvent) and therefore 
the sum of the molar fractions is approximately: 

X *,= I >,[ 1 + (*',- 1 )«,] [ 5 - 31 ] 

s s 


For each dissociable solute, we see the emergence of a factor i s defined 
by: 


i s =[\+(v s -\)a s ] 


[5.32] 


This factor i s is known as the Van ’t Hoff factor. 

The relation found by ebullioscopy therefore involves the mean activity 
coefficient. According to expression [5.26], it is written thus: 


In Y± = 


A v h 0 AT( Eb ) 

n T 2 



[5.33] 


If the solution contains only one electrolyte, then relation [5.33] can be 
simplified to: 


In r± 


A v h 0 Af E b) 

R T 2 


-ii 


(L) 


[5.33a] 


Similarly, on the basis of cryoscopic measurements, relation [5.28] gives 
us: 


In r± 


Apho ^T^f) 

R T 2 



s 


[5.34] 
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In the case of the dissolution of a single electrolyte, relation [5.34] can be 
simplified to: 


In y± = 


A F h° 0 AT (F) 

R T 2 

S ^ 1 0(F) 


-i.x. 


{L) 


[5.34a] 


Thus, we have two equations [5.33a] and [5.34a] and two unknowns: the 
coefficient of dissociation a s and the mean activity coefficient y ± . 
Obviously, this coefficient will be obtained for a mean temperature between 
the solvent’s freezing and boiling points. 

NOTE - In the case of strong electrolytes (a s = 1), the Van ’t Hoff factor is 
equal to the number of ions v s supplied by the dissociation of that 
electrolyte. 


5.3.3. Use of the measurement of osmotic pressure 

Remember that the osmotic pressure is the pressure difference between a 
solution and its solvent, separated by a semi-permeable membrane, which 
allows the solvent, but not the solute, to pass through (see Figure 5.3). 


P 

4 



Figure 5.3. Osmotic pressure 


Taking account of the variation of the chemical potential with the 
pressure, we can show that for sufficiently dilute, wherein the partial molar 
volume of the solvent in the solution can be treated as the same as its molar 
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volume in the pure state, the equality of the solvent’s chemical potential in 
both phases gives us the relation: 


In 7o = 


17 V 0 
R T 


[5.35] 


In this relation, 77 is the osmotic pressure, i.e. the difference between the 
pressure P over the solution and the pressure P tl over the solvent (see 
Figure 5.3). 

If the solution contains several solutes with the molar fraction x s , this 
expression can also be written as follows, for a dilute solution: 

ln r 0 =-fp L +Z^ C 5 - 36 ] 


In order for the accuracy to be sufficient, the osmotic pressure - i.e. the 
height h measured in Figure 5.3 - must be significant. This result is easy to 
obtain with solutions of polymers, for which numerous semi-permeable 
membranes are available. It is more difficult to find a good membrane for 
smaller molecules of solutes. Nevertheless, often, the measurement is 
significant because, for instance, for an aqueous solution of sodium chloride 
of 0.1 molar concentration, at 25°C, the osmotic pressure reaches 
4.61 x 10 5 Pa. 


5.4. Measuring the activity on the basis of solubility measurements 

The solubility which results from an equilibrium of a component between 
two phases - one the pure phase of one of the components, and the other the 
solution - can also be used to determine the activity or the activity 
coefficient of a component in a solution. In general, it pertains mainly to 
solutes, but it is also applicable to solutions for which the concept of a 
solvent is no longer clear. The two phases including the component in 
questions are usually the pure solid phase on the one hand, and the liquid 
solution on the other. 

Nevertheless, we shall distinguish between two cases, depending on 
whether the solution is molecular or ionic, i.e. whether the species in 



174 Modeling of Liquid Phases 


solution is identical to that of the solid, or whether it has undergone 
alterations as it dissolved. 


5.4.1. Measuring the solubilities in molecular solutions 

The molecule of the component in pure phase and that of the solution are 
identical; thus, there is a simple equilibrium between the two phases, so the 
chemical potentials of the component concerned are equal for both phases. 
Thus, we find the following relation, which is valid regardless of the 
convention chosen for the solution: 


In fi 


hl) _ A F h, 


R 


V T i(n 


T 


■ ln d-2>f) 




[5.37] 


If the solution has a very high content of component 0 (as is the case with 
a solvent), this relation is simplified, for that solvent, to: 


In 

R 


1 


T 

V 0(F) 


T 


+ 


Z- 

j* 0 


.(£) 


[5.38] 


Thus, relations [5.37] or [5.38] enable us to calculate the activity 
coefficient of a component if we know the composition of the solution, the 
temperature of solubility equilibrium chosen, the melting point of the pure 
component in its natural phase and its standard enthalpy of melting at the 
temperature in question. 


5.4.2. Measuring the solubilities in ionic solutions 

We shall only envisage the case of complete dissociation of the 
electrolyte, which is true for all salts and most other species with low 
solubility at the solubility equilibrium, because then the solutions are very 
dilute. 

The species in solution (the ions) result from this dissociation upon 
dissolution. The equilibrium between the solid and the solution, therefore, is 
characterized by a solubility product K s . If S denotes the solubility measured 
by the amount of solid which has passed into solution (as opposed to the 
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amount of ions in the solution), and if the dissociation of a molecule of solid 
yields V + cations and V_ anions, we can show that writing the equilibrium 
gives us the following relation: 


/ ~ionst . 

S = - [5.39] 

r± 


y ± is the mean activity coefficient and the constant, which depends only 
on the temperature, is linked to the solubility product by: 


^jonst. 


K, 


\ 


v V *v K 

V K + - / 


l/(v++v + 


[5.40] 


If we switch to a logarithmic system, relation [5.39] is rewritten as: 

In y ± = In C°” s ' - In S [5.41] 

The mean activity coefficient tends toward 1 when the ionic strength of 
the solution tends toward 0. Thus, we deduce relation [5.42]: 


lnC°" s ' =hm(lnS) 


[5.42] 



Figure 5.4. Solubility of silver nitrate as a function 
of the ionic strength (reproduced from [POP 30]) 
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Thus, by measuring the solubility S at different ionic strengths / and 
extrapolating the curve at an ionic strength of zero (Figure 5.4), we obtain 
the value of the constant C’" sL in relation [5.41], because then the mean 
activity coefficient tends toward 1. From this we deduce the value of the 
solubility product K s by relation [5.40], and that of the mean activity 
coefficient, regardless of the value of the ionic strength, by relation [5.41]. 


5.5. Measuring the activity by measuring the distribution of a solute 
between two immiscible solvents 


If we know the activity coefficient y^ a 1 of a component in a phase (a), for 
example, and the partition coefficient K.\ a/I) of that component between that 

phase a and another phase /? which is immiscible with it, then it is easy, by 
describing the equality of the chemical potentials of the component i 
between the two phases, to obtain the other activity coefficient y'/ 1 ’ using 
the relation: 


'> v </?) 


[5.43] 


Thus, it is only necessary to know the molar fractions x. 0) and x'/ l] of 
the component in the two phases at equilibrium and of the activity 
coefficient in phase a to obtain . Relation [5.43] is independent of the 
reference state chosen. 


5.6. Activity in a conductive solution 

If the solution at hand is able to conduct electrical current, the panoply of 
methods that can be used is enriched with a set of electrochemical methods. 
These methods can be used in the case of ionic solutions used as electrolytes. 


5.6.1. Measuring the activity in a strong electrolyte 

In a strong electrolyte, we may be led to measure either the absolute 
activity of an ion or the mean activity of the electrolyte. 
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5.6. 1 . 1 . Measuring the absolute activity of an ion 

The measurement of the absolute activity of an ion of a strong electrolyte 
(salts, strong acids and bases) is based on the complete dissociation of that 
electrolyte into ions, and uses the measurement of the electromotive force of 
a cell involving that ion. 

In principle, we use the fact that electrode potentials are linear functions 
of the logarithm of the ion activity if the cell’s behavior is reversible. 

Thus, it is sufficient to construct a battery formed of an active electrode 
of the ion under study - the measuring electrode - and a reference electrode. 
In general, the liquids in which the electrodes are bathed are connected by a 
siphon, filled with a concentrated electrolyte in order to minimize the 
junction potentials. 

For example, in order to measure the activity of the hydrogen ion, we can 
construct a system containing a mercurous-chloride electrode (the reference 
electrode) and a hydrogen electrode (the measuring electrode), connected by 
a junction siphon (see Figure 5.6). 

The potential, in this case, is measured at 25°C, and obeys the relation: 

£ abs =0.252 -0.06 In |H + 1 [5.44] 



Thus, by measuring this potential, we are able to measure the activity of 
the solutions in terms of protons. 


178 Modeling of Liquid Phases 


A variant of this method is to create a concentration cell by taking two 
identical electrodes, immersed in two solutions with different activities, one 
of which is known. 

For example, we immerse two silver electrodes in two solutions of Ag + 
ions - one in which the activity of the silver is known (solution 2) and the 
other being the solution in which we wish to measure the activity in silver 
ions (solution 1). The potential of the cell thus constituted will be: 

£abs =e% +0.°61n|^g + | 2 -e% - 0.06 In [5.45] 
Thus: 


= 0.06 In 



[5.46] 


This relation enables us to measure the activity of the unknown solution 
without needing to know the value of the standard potential of silver e l ’ v . 

In practice, the true value of the activity of an ion is difficult to find 
because of the error introduced by the junction siphon. Although the 
assembly attempts to minimize the junction potential, our calculation 
supposes that the potential is null, and the true value of the activity 
coefficient of the ion is tainted with an error. 

5.6. 1.2. Measurement of the mean activity coefficient of a strong electrolyte 

In order to obtain more accurate results, it is preferable to use cells with 
no junction, but in that case, the method of measuring the activity 
coefficients of strong electrolytes will yield only the mean activity 
coefficient. 

If we know the standard potential of a carefully-chosen cell (which may 
be determined experimentally or calculated on the basis of the electrode 
standard potential tables), then we can calculate the mean activity coefficient 
at any concentration. 
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We shall demonstrate the method with the example of measuring the 
mean activity coefficient of a zinc chloride solution of concentration c. For 
that purpose, we create the following cell with no junction: 

Zn / ZnCl 2 (solution) / AgCl (dissolved), Ag 

The cell’s emf would be: 


E^=E°-— ln|Zn ++ ||cr| 2 [5.47] 

This involves only the product of the activities Zn ++ Cl“ rather than 
the individual activities of the ions, which would be: 

|Zn ++ 1 = [Zn] y Zn „ = c/ Zn++ and \ci~ | = [_Cl~ ] y a _ = 2cy cr [5.48] 

By combining relations [5.47] and [5.48], the emf of the chain becomes: 

!»■ 4c’ - ^ In r Zn .f cr [5.49] 

Flowever, the mean activity coefficient is defined by: 

rl = r^fa- [ 5 - 5 °] 

Relation [5.49] can thus be written in the form: 

E ^= E °- — ln4c 3 -— lny + [5.51] 

abs 2F 2F 

If we know the standard cell potential E° and the concentration c, we can 
deduce the mean activity of the ions using the relation derived from the one 
above: 


[5.52] 
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Let us stress the fact that it is advisable to create cells without junctions, 
because no matter how perfect those junctions may be, they cause significant 
errors of up to 50% on the value of the cell emf. 


5.6.2. Determination of the mean activity of a weak electrolyte on the basis 
of the dissociation equilibrium 

Suppose that we know the dissociation constant of the weak electrolyte 
AH and its dissociation coefficient a at the chosen concentration c. The 
dissociation constant is written: 


K _ etc Y± 
l — OC y AH 


[5.53] 


Hypothesize that we are dealing with a sufficiently dilute solution of 
electrolyte AH. We know that when the concentration decreases, solutions of 
neutral molecules approach the ideal state much more quickly than do those 
containing ions. Hence, we shall set y AU = 1 . Then, we immediately obtain: 


iK(l-a') 

X + =J , ’ [5.54] 

V ore 

This simplification means that we now only have one unknown value - 
the mean activity coefficient of the ions in the solution. Such an 
approximation is correct provided the concentration of electrolyte is no 
greater than 1 O' 1 - 1 0 2 moles per liter. 
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Appendix 1 


Statistical Methods of 
Numerical Simulation 


The purpose of this appendix is not to give readers a grounding in 
statistical simulation, but rather to enable future thermodynamicists to 
conduct a dialog with simulation specialists. We wish to demonstrate the 
physics which lies behind each simulation, and its reach, so that it is possible 
to appreciate the thermodynamic results of a simulation that the reader has 
conducted in complete collaboration with the specialists. 

The objective of a simulation is to generate particle motion by using 
appropriate algorithms and to obtain an adequate distribution function, and 
thus, the macroscopic properties. In thermodynamics, we use this method to 
calculate the energy of interaction of a collection of molecules, the part of 
configuration of the translational partition function and the radial distribution 
function. 


A.1.1. The physical bases of simulation 

A simulation, whatever the method employed, is based on physical 
hypotheses. In the case of the molecular simulation, there are two 
fundamental hypotheses: 

- the first pertains to the link between the internal energy of the set of 
molecules and the interactions which exist between them. One of the most 
commonly chosen hypotheses is to consider only the sum of the pairwise 
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interactions between the molecules, excluding interactions involving more 
than two substances. In the case of the internal energy of configuration, it 
obeys the relation: 


U I {l,2,...N) = f j £ iJ (r iJ ) 

i<j 


[A.1.1] 


Thus, we can calculate the configuration part of the canonical partition 
function of translation, which is: 


h =Zexp 


*</ 


<{ r u) 


k B T 


[A. 1.2] 


We can also obtain the radial distribution function given by: 

1 N N 


s( r ) = -^YL S iA r - r i,j) 

ply i = i j = i 


[A. 1.3] 


-the second fundamental hypothesis is the choice of the potential 
function of interactions between molecules. In the hypothesis chosen in 
relation [A. 1 . 1 ], it is the choice of the function e t j ( r. ; ) . 

A variety of laws have been encountered for the potential of attraction: 

-the Lennard-Jones potential in \/<f, which relates to the van der Waals 
forces between molecules and stems from electrostatic forces between 
electrical dipoles. This potential is used for real gases and liquids; 

- the Coulombian potential in 1 Id, stemming from the electrostatic forces 
between two points or two conductive spheres, used between ions in an ionic 
solution; 

- the ion-dipole potential in l Id 2 of electrostatic nature between a 
conductive sphere and an electric dipole, which is useful, for instance, for 
examining the solvation of ions in a dipolar medium such as water; 

- the potential for interaction between a molecule and a surface, in 1 Id 2 , 
used for physical adsoiption, which is the resultant of the different Lennard- 
Jones potentials. 
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For the repulsion potential, we normally use a potential in d xln , frequently 
with n =12, which expresses the repulsion to the interpenetration of the 
electron clouds, which is therefore used for ions as well as molecules. To 
avoid an integral which diverges as the distance tends toward zero, we often 
adopt the approximation of a hard sphere, which gives an infinite potential 
for a protection distance less than or equal to D which is the minimum 
approach distance. 

Of course, the overall potential requires the combination of potential of 
attraction and a potential of repulsion in order to represent the possibility of 
a position of equilibrium. 


A.1.2. Construction of the sample 

If we consider a system containing, for example a mole of liquid, this 
represents around 10 24 molecules, which represents 10 24 x 10 24 /2, which is 
around 10 48 couples of molecules that need to be taken into account in 
relation [A. 1.1], for all combinations of relative positions of those couples 
(the configurations of the system). There is no computer system, no modem 
memory, that is capable of storing such a phenomenal amount of 
information. Therefore, we need to perform the simulation on a much 
smaller sample, whilst remaining faithful to the behavior of the real system. 
The sample adopted, therefore, needs to satisfy two fairly contradictory 
conditions: 

- it should not contain too many elements in order to be manageable; 

- it should contain enough elements for it to still be representative of the 
system. 

In order to resolve this dilemma, we use a variety of techniques: 

- limiting the distance over which the interactions between molecules are 
taken into account, by truncating the potential function; 

- eliminating the edge effects, which can become very significant, 
relatively speaking, in a small system; 

- testing the validity of the calculation at every step to limit its duration. 
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A.l.2.1. Truncation of the potential function 

The curve showing the potential function generally tends toward zero 
when the distance between two elements increases, as is shown by 
Figure A. 1 . 1 . This means that, if we accept a slight degree of error, we can 
take account only of a molecule’s interactions with its fairly close neighbors, 
that is, if the distance r . between two elements is less than a given distance 

r c , called the cutoff distance (Figure A. 1.2), so: 
if r i,j ^ r c then £ = 


[A. 1.4] 


if r , > r then £ = 0 

We then compile what are called Verlet lists which, for each element in 
the sample, give the list of all the elements j which are at a distance less 
than a certain minimum distance r m (Figure A. 1.2), so as to be able to 
perform multiple steps of calculation without changing the list relative to a 
particle. 



Figure A. 1.1. Lennard-Jones potential and cutoff radius 


The minimum radius r m is such that: 



[A. 1.5] 


In this formula, k is the number of times we use the list without changing 
it, At is the step increment (see section A. 1.3.1). It is advisable to take the 


Appendix 1 187 


size of the sample for the calculation as that which contains N smp i elements, 
which is given by: 

Ajr 

K^=~prl [A. 1.6] 

Here, p is the density of the real system in terms of molecules. 



A.l.2.2. Limitation of edge effects 

When we choose a set with a small number of elements, the edge effects 
are very significant, meaning that a high proportion of the elements are near 
to the boundary of the system, and this poses two problems: 

- during the course of a displacement, an element may stray out of the 
domain, and therefore we are no longer working with a constant number of 
particles. We remedy this problem by using what is known as the periodic 
boundaiy condition ; 

- an element near to the edge of the domain no longer inhabits its 
“normal” environment. We correct this effect by using the minimum-image 
convention. 

A. 1 .2.2. 1 . Periodic boundaiy condition 

Let us consider that the system in question is delimited by a cube of side 
length a. We can construct the cubes adjacent to that initial cell by exactly 
reproducing the composition of the initial cell. Figure A. 1.3 offers a 
2D representation of the stack thus created. The three Cartesian axes can be 
used to write, for example, that the central cell in our study is such that its 
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boundaries have the abscissa values X=al2 and X=-a/2, and similarly in 
the other two directions. During the displacement of the elements, one of 
them (A) near to the boundary may leave the central cell, crossing that 
boundary. So as to work with a constant number of elements, we shall use 
the Born and von Karman method, which consists of bringing in a 
corresponding element (B) across the opposite boundary, i.e. setting the 
following condition, for example, along the abscissa axis: 

if X: <—a/2 then x=x+a 

1 11 [A. 1.7] 

if Xj > a / 2 then x ] =x--a 


• c 4r • cf 


O' 


.ir 


• o 


-►.v 


• O' 


-• • 

• Ol 




o' 


d e • o*" • c*'' 




Figure A.1.3. Born-von Karman periodic boundary condition 


A. 1.2. 2. 2. Minimum-image convention 

If an element A is too close to the border of its environment, there is a 
risk that it will not be complete. In order to correct this effect, we apply 
the minimum-image convention, which consists of including in the 
Verlet list for that element A, and in the calculation (Figure A. 1.4), the 
elements neighboring A, situated in adjacent copies of the calculation cell, 
such as B. 
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Figure A.1.4. Minimum-image convention 


A.l.2.3. Estimation of the duration of the calculation 

All of these calculation methods work on a step-by-step basis. A step has 
an increment, denoted by At. In general, we begin with an initial 
configuration - e.g. the compact stack represented in two dimensions in 
Figure A. 1.5 - and then, by successive displacements, at each step, we 
obtain a succession of configurations. One problem that arises during 
the calculation is knowing the number of steps after which the calculation 
should be halted in order to be representative of the state of the true system. 



Figure A.1.5. The initial compact hexagonal stack 


To determine this number of steps, we use a coefficient called the self- 
correlation coefficient on g(r), for example. 

The self-correlation coefficient measures how the value of x at step t+ At 
is connected to the previous value of x at step t; or, if readers prefer to think 
of it this way, it measures the influence of the previous state on the present 
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state. This function is calculated by relation [A. 1.8], which gives a curve that 
begins at 1 in the initial state and tends toward zero after a relatively long 
period of time (see Figure A. 1.6). We need only set the acceptable error to 
determine the stopping point of the calculation. 

M 

Y J [x{t k )-(x)][x(t k+ ^)-(x)~\ 

C(t) = - * [A- 1-8] 

Z [*(**)-(*)] 2 


cm 



A.1.3. The main calculation methods 

There are two types of calculation methods: static methods, which lead to 
a state deemed to be stable, and dynamic methods, which show the system’s 
evolution over time until it reaches a stable configuration. 


A.1.3. 1. The Monte-Carlo method 

The Monte-Carlo method is a static statistical simulation method. 

In this method, we begin with an initial configuration and the total energy 
of the system U\ is calculated on the basis of relation [A.1.1]. Then, each 
element is shifted in the three directions of small amplitudes Ax, Ay and Az. 
This amplitude is chosen as around 10% of the molecular diameter. Thus, 
each element has a new position, which gives a new configuration, for which 
we calculate the new energy U 2 . If U 2 is less than U\, the new configuration 
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is accepted and the coordinates of the elements are saved to memory. If the 
new energy U 2 is greater than U\, we calculate the ratio fl defined by: 


f 

n = exp 

V 


U^_ 
k J 


f 

/exp 

V 


Uy_ 

KTj 


[A. 1.9] 


This ratio is compared with a random number X between zero and one. If 
X is smaller than II, the new configuration is refused and the first 
configuration is saved to memory again. If X is larger than n, the new 
configuration is accepted and saved to memory. Then, we carry out a new 
calculation step. All the elements are again shifted by the same distances as 
before, and the new value of n is compared to a new random number X’, and 
so on. Thus, we find around 10 5 to 10 6 configurations. 

The increments, Ax, Ay and A z, of displacement must not be too large, 
otherwise few solutions are accepted, nor too small, which would lead to too 
high a number of acceptable solutions. 


A.l.3.2. The molecular dynamics method 

In this method, the evolution over time of the system with of N elements 
is studied on the basis of the fundamental law of dynamics: 


F — mid [A. 1.10] 

The calculation step, in this method, is a period of time. The forces on 
each element are calculated at each step, with: 



Z 

. 1< J 


[A.l. 11] 


For this, we use an appropriate algorithm, of which one of the best known 
is attributable to Verlet. In this algorithm, the displacement of a particle r at 
time t + A t is chosen on the basis of the previous two displacements at time t 
and t — At. Thus, we obtain the positions and velocities of all the elements at 
each step of the calculation. Based on those data inscribed in memory, we 
can calculate the radial distribution function using relation [A. 1.3]. 




Appendix 2 

Reminders of the Properties of Solutions 


A.2.1. Values attached to solutions 

Consider a solution comprising N components which contains n t moles of 
component 1 , n 2 moles of component 2, n t moles of component i, etc. 

The molar fraction of component i is x„ defined by: 

N 

= y ‘ where 0 < x t < I and Vx, =l [A.2.1] 

& 1=1 
i = 1 

We use /jj to denote the chemical potential of component i. 

This solution is completely characterized by the knowledge of the 
function the Gibbs energy G, which is written as follows, in differential 
form: 

N 

dG = -SdT + PdV + Y J M i dn i [A.2.2a] 

1=1 

In order to find this function, we only need to know the chemical 
potentials, because they are the partial molar Gibbs energies, and we have: 

G = Y,m [A.2.2b] 
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A.2.2. Peculiar values and mixing values 

Generally speaking, an extensive value - e.g. the volume of a solution - 
is not the simple sum (weighted or otherwise) of the volumes of substances 
introduced when making up the solution. Usually, an extra term is added 
which is attributable solely to the mixing of the substances, called the 
corresponding mixing value. 


A.2.2. 1. Definitions 

Consider a partial molar value Jj . In a uniform solution, that value is a 
function of the variables temperature and composition. If we examine a mole 
of the overall solution, taking as composition variables the molar fractions x h 
it is always possible to decompose that function into the sum of two 
functions such that one depends only on the temperature (we ignore the 
influence of pressure on a liquid solution) - this is the peculiar function 

pec 

Jj - and the other is a function of all the variables of temperature and 

— mix 

composition, and is known as the mixing function ./, . This decomposition 

is written thus: 


— — pec — mix 

J i =J i ( T ) + J i (T,X\X 2 ---X N ) [A.2.3] 

As a peculiar value, we choose the molar value of the pure substance i in 
the same state of agglomeration as the solution. The decomposition [A.2.3] 
is then written: 


mix 

J, =f(T) + J, (T, Xl x 2 ...x N ) [A.2.4] 

ji denotes the molar value of the pure property: 

The operation [A.2.4] can also be applied to any linear combination of 
partial molar values, and in particular to the combination: 

N 

Ji=Yj n i J i 

/=i 


[A.2.5] 
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Mixing 

value 

Definition 

Expression 

Gibbs 

energy 

G'm X = IaU - £/°) [A.2.6] 

1=1 

N 

G™ x = RT J] Xj In Xj [ A.2.1 ] 
1=1 

Helmholtz 

energy 

HT = - h?) [A.2.8] 

i=l 

= 0 [A.2.9] 

Entropy of 
mixing 

S™ x = Z^(^-«1°) [A.2.10] 

i=l 

Sm* =-RT J X i lwc i [A.2.1 1] 

i=l 

Mixing 

volume 

C“ = Zv(h-V,°) [A. 2.12] 

1 = l 

C“' = 0 [A.2.13] 

Specific 

heat 

capacity at 

constant 

pressure 

Cp"n X = Za(Cp ; -CP')[A2.14] 
1=1 

N 

Cp = X4(1) [A.2.1 5] 

i=l 


Table A.2.1. Definitions of the mixing properties and values for a perfect solution 


We use to denote the molar value of mixing of the solution, defined 
by: 


j: 


N 

= 1 * 

i= 1 




[A.2.1 6] 


The second column in Table A.2.1 gives the definitions of several molar 
values of mixing. 
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A.2.3. Characterization of the imperfection of a real solution 

A perfect solution is defined as a solution in which the chemical potential 
of each of the components obeys the relation: 

Mi — Si + RTlnx,- [A 2 17] 

From this definition, we deduce a certain number of classic properties of 
the perfect solution. Table A.2.1 gives a few molar values of mixing for a 
perfect solution. 

Two methods have been distinguished. What they have in common is that 
they characterize a real solution by its difference from a perfect solution: 

- the first method, from Lewis’ school of thought, introduces the activity 
coefficients; 

-the second method uses the excess values, in particular the excess 
Gibbs energy. 


A.2.4. Activity coefficients 

Lewis refers to the expression [A.2.1 7] of the chemical potential, 
attempting to preserve its form. In order to do so, he introduces an activity 
coefficient % which is a function of the temperature and of the composition 
of the solution, writing the chemical potential of a component i of the 
solution in the form: 

/U- = jl] + RTIn y i x i [A.2.1 8] 


A.2.5. Activity coefficients and reference states 

The product of the activity coefficient and the molar fraction is called the 
activity of the component i in the solution. The activity in a real solution 
plays the same thermodynamic role as the molar fraction does in a perfect 
solution: 


a, = /jXj 


[A.2.1 9] 
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The reference chosen for convention (I) is components in the pure state in 
the same state of segregation as the solution (this reference is called the 

pure-substance reference), in which case, y\ l ] = 1 . We chose to use it in 

section A.2.2 for the mixing values. In these conditions, the chemical 
potential of the reference state is the molar Gibbs energy of the pure i, and 
then the chemical potential is written: 

Mi = g? {T)+RT In f ) x i [A.2.20] 

This reference is mainly used when all the components in a solution play 
the same role, and in particular, have comparable molar fractions. For 
example, this convention is chosen when the data cover a broad spectrum of 
composition, possibly ranging from one pure substance to another. 

Convention (II), called the infinitely-dilute solution reference, 
distinguishes, among the components of the solution, that (or those) present 
in a high proportion, called the solvent(s), and those which are present in 
lower proportions, called the solutes. The reference state is different for 
these two categories of components: 

- for a solvent, we choose its pure state (in the same state of segregation 
as the solution) as a reference, and therefore its chemical potential will 
obey relation [A.2.20]. For a solvent, convention (II) is the same as 
convention (I); 

- for a solute, the reference state is an imaginary solution in which all the 
solutes are infinitely dilute. The reference chemical potential is therefore that 

of the solute in this imaginary solution, and is written as juf for a solute 5. 
The activity coefficient is equal to 1 in this imaginary solution and the 
chemical potential will then be written, for all solutes (accepting that the 
activity tends toward 1 if the molar fraction tends toward the composition of 
the reference state): 

M s ~ M7 ( T) + RT In '/ s in x s [A.2.2 1] 

The activity coefficients of a solute in reference (I) - pure substance - 
and (II) - infinitely-dilute solution - are linked to one another because the 
chemical potential of the solute does not depend on the convention chosen. 
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From this, we deduce: 


ff _ M7(T)-g°s(T) 

y/ n R T 


= ln K ih 


[A.2.22] 


Thus, the constant K iH linking the activity coefficients expressed in the 
two conventions: the pure-substance reference and the infmitely-dilute 
solution reference do not depend on the composition of the solution, but 
instead depends on the temperature by means of (amongst others) the 
chemical potentials of the reference states. The value K iH is called Henry’s 
constant. This constant can be determined, in particular, on the basis of the 
measured vapor pressure at equilibrium between the solution and the vapor 
phase (see section 5.2.2). 


A.2.5.1. Relation between the activity coefficients of the components of a 
solution 

Let us place ourselves in the context of any given convention. The 
Gibbs-Duhem relation applies to chemical potentials which, at constant 
temperature, obey the relation: 

£jc,d/4=0 [A.2.23] 

1=1 

From this, we deduce: 

^ x t d In Yi = 0 [A. 2. 24] 

/=! 


Thus, if, for each composition of the solution, we know the activity 
coefficients of all the components in a convention, except for one of them, 
relation [A.2.24] can be used to calculate the unknown activity coefficient in 
the same convention (see section 5.1.1). 
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A.2.5.2. Influence of the different variables on the activity coefficients 

A.2.5.2.1. Temperature 

If H{ is the partial molar enthalpy of component i in the reference state 
at temperature T, we can show that: 


d In y i _ //" - //. 
dT ~ R T 2 


[A.2.25] 


This relation [A.2.24] gives us the variation of the activity coefficient at 
constant composition when the temperature changes. 

The numerator in expression [A.2.25] reveals a difference which 
represents the enthalpy of transport, at constant temperature, of a mole of i, 
from the real solution under study to the solution in the reference state. We 
must imagine that the transport takes place with a large amount of solutions 
so it does not alter the compositions of the two solutions. If we adopt 
convention I, it is clear that this difference is the opposite of the enthalpy of 

mixing (because then = hf ) . 

A.2. 5.2.2. Influence of the composition on the activity coefficients 
If a solution is perfect, it obeys Raoulf s law, which is written as: 

a\ ,) =x i [A.2.26] 

This law is represented by Raoulf s straight line, which is the first 
bisector in the system of axes ( a^\x t ). 


An ideal dilute solution is defined by an activity coefficient in reference 
(II) equal to 1 . Its activity coefficient in reference (I) will therefore be given 
by relation [A.2.22]. Thus, the activity of i in convention (I) in an infinitely- 
dilute ideal solution is: 


af = K in x, 


[A.2.27] 
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This is Henry’s law. This law is represented by Henry’s straight line, 
whose slope is K m in the system of axes ( ,x t ). Henry’s constant K iH is 

independent of the composition and has a very clear value when the 
temperature is constant (see relation [A.2.22]). 




Figure A.2.1. Activity curve of a component as a function of its molar fraction for 
a) a component with positive deviation and b) a component with negative deviation 


For a real solution, the activity, in the pure-substance reference, of a 
component tends toward 1 when its molar fraction tends toward 1 . It tends 
toward 0, along Henry’s straight line, when its molar fraction tends toward 0. 
Thus, the curve showing its activity, in the pure-substance reference, as a 
function of its molar fraction, is tangential to Raoult’s straight line at the 
point x = 0. It is tangential to Henry’s straight line at the point x = 1 . 

A solution is said to show positive deviation if its activity coefficient in 
the pure-substance reference is greater than l(y () >1), the curve is above 
Raoult’s straight line. This solution is said to show negative deviation if its 
activity coefficient in the pure-substance reference is less than l(y (I) <1), 
with the curve below Raoult’s straight line. 

Besides, we can show the inter-relation between the activity coefficients 
of two components in a solution: 

3 In y t _ 3 In y k 
dx k 3x, 


[A.2.28] 
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This relation stems from a similar relation between the chemical 
potentials, which is the consequence of the symmetry of the characteristic 
matrix. 


A.2.6. Excess values 

We can show that the activity coefficient of a component in a perfect 
solution in the infmitely-dilute solution reference is also 1 at all temperatures 
and in all compositions. 

tfo = ^) = i [A.2.29] 

NOTE.- Henry’s constant for a perfect solution has a value of 1. 

Consider an extensive property J of a solution and use J ^ to denote the 
value of that property in given conditions of temperature and pressure if that 
solution were perfect. We speak of the excess value of J, denoted by J xs 
and defined by: 

J xs =J-J P f [A.2.30] 


J xs is indeed characteristic of the difference between our real solution 
and a perfect solution. 


As this value is extensive, it has corresponding partial molar values 
pertaining to each component in the solution, which are defined by: 


J xs = 


dr 

v dn > 


[A.2.31] 


Table A.2.2 presents a few excess values, along with the corresponding 
excess partial molar values. 

This table also demonstrates the equivalence between the two approaches 
to the modeling of real solutions: Lewis’ approach based on the activity 
coefficients, and the approach based on the excess values. 
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Property 

Excess molar value 

Excess partial molar value 

Gibbs 

energy 

G™ =Rrf> f ln y\ I] [A.2.32] 
i = 1 

Gf = R T In 

[A.2.33] 

Entropy 

n r 

S” = R Z x i 

i=i L l 

din y\ I) 

T ' 1 in V; 

dT 1 

J _ 

[A.2.34] 

Sf s = R dT 

\r ln rt I] j 

[A.2.35 

1 

Enthalpy 

N 

ttXS _ p V'' 

i=l 

dT 

[A.2.36] 

r? o d In y< 7) 

H xs = rt 2 h 

dT 

[A.2.37] 

Specific 

heat 

capacity 

N 

c x ; = ri x 

i=l 

dT 2 

0 dln r\ I] 

l dT J 

[A. 

2.38] 



Table A.2.2. Excess values and excess partial molar values 


A.2.7. Ionic solutions 

Ionic solutions exhibit several peculiarities in comparison to molecular 
solutions. The main source of these differences is to be found in the 
existence of inter-ion forces that are much stronger than those that exist 
between molecules - particularly over a long distance. 


A.2.7.1. Composition of an ionic solution and reference state 

To express the composition of the ionic solutions, we usually use the 
concentration (or molarity) c t of a component i in a phase. This is the 

quotient of the amount of that component by the volume of the phase, so: 
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c, = 


>h_ 

V 


[A.2.39] 


The concentration is often expressed in moles/liter. For solutions in 
which the proportion of one of the components is much greater than the 
others (often that component is called the solvent), we frequently treat the 
volume of the phase as the same thing as the volume of the component 
introduced in the highest quantity when making up the solution. Its 
concentration is often considered constant, with the addition of a small 
amount of another component. 

The reference state is usually defined by convention (III), which makes 
the distinction between solvent and solute: 

- for the solvent, the reference (III) convention is identical to 
reference (I) - the pure substance - and therefore the chemical potential of 
the solvent is always given by relation [A.2.20] ; 

- for the solute, we shall agree to choose that the activity coefficient of an 
ion tends toward 1 as its concentration tends toward zero, which gives us the 
following convention: 

y s ->l if c, -» 0 [A.2.40] 


A.2.7.2. Chemical potential of an ion 

In spite of the difficulty, in the case of ions, of defining a derivation in 
relation to one of the components whilst preserving the values of the other 
constants, we can show that it is exactly as if we had chosen the usual 
relation to express the chemical potential: 

M f ! + R/ ln YjCj [A.2.41] 

The potential //,° is, of course, that of the reference state. 
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A.2.7.3. Relation between the activities of the ions and the overall activity 
of the solutes 

Consider the dissolution of an electrolyte with the chemical formula 
A,, B v _. Suppose that it is a strong electrolyte; it dissociates completely. Its 
Gibbs energy is written as: 

G = Oq/Jq + n s + v_//o j [A. 2. 42] 

However, we know that the chemical potential of the solute is also: 

Ms=v+Ma+V-Mb [A.2.43] 

From this, we deduce: 

r s c s ={rAC A ) v+ {r^ B ) v - [a.2.44] 

This equation links the molecular point of view about the solution to the 
ionic point of view. 


A.2.7.4. Mean concentration and mean activity coefficient of the ions 


The methods for precisely measuring the activity coefficients (see 
section 5. 6. 1.1) are incapable of giving us the activity coefficients of the 
individual ions. Thus, it has proved helpful, for an electrolyte A v _, to 

introduce the concept of the mean activity coefficient defined by: 


r± = 


= ir+‘r- v -] 


'{v + +v_ 


[A.2.45] 


The mean activity coefficient obeys the same convention as the 
individual activity coefficients. 

We can also define a mean concentration by a similar relation. If c is the 
molar concentration of solute, the concentrations of the different ions (for 
fully-dissociated strong electrolytes) would be: 


c + = v + c and c_ = v_c 


[A.2.46] 
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and the mean concentration would be: 




l/(v + +V_)) 


[A.2.47] 


A.2.7.5. Activity coefficient of an individual ion 

We have already mentioned how difficult it is to determine the individual 
activity coefficients for the ions by experimentation. It is obvious that if we 
were able to determine the activity coefficient of a single individual ion, then 
little by little we could - on the basis of the mean activity coefficient, which 
we are able to measure (see section 5.6. 1.2) - deduce that of all the other 
ions. In order to do so, Maclnnes exploited the fact that in potassium 
chloride, the chloride ion and the potassium ion have the same charge in 
absolute value, have the same electron structure and therefore the same size, 
roughly the same mobility and masses that are not hugely different. 
Maclnnes deduces from this that they should have approximately the same 
activity coefficient. In light of that result, we deduce, from the mean activity 
coefficient of potassium chloride, that of the separate ions, using the relation: 

Y cr = r K+ = r± [A.2.48] 

By the same hypothesis, the results obtained for a set of ions show that 
monovalent ions have approximately the same activity coefficient, and that 
the activity coefficients are primarily influenced by the electrovalences of 
the ions and by the presence of other ions, which gives rise to the concept of 
ionic strength. 


A.2.7.6. Concept of ionic strength 

The activity coefficient of the ions is influenced by the presence of other 
ions. Experience shows us that the important value is the ionic strength, 
defined by: 

I= \H c i z i t A - 2 - 49 ] 


The sum is extended to all the ions present in the solution. 
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Reminders on Statistical 
Thermodynamics 


The purpose of this appendix is to recap, but not demonstrate, a few 
results of statistical thermodynamics which are used in this second volume. 
The detail and expansion of these concepts are presented in the first volume 
of the collection [SOU 15]. 

We know that a microscopic approach to a phase considers it as a 
collection of molecules-objects whose energies are distributed in accordance 
with a statistical law. The state of a collection of molecules-objects is 
constantly changing but, over time, the collection goes through a certain type 
of distribution in which the molecules-objects may be in a variety of 
different states. We use the term number of complexions for the number of 
distributions of molecules-objects between the different states that they are 
liable to assume. Out of all the possible distributions, there is one that 
corresponds to the maximum number of complexions. Boltzmann’s principle 
accepts that the number of complexions corresponding to the most probable 
type of distribution is practically equal to the total number of complexions, 
and vice versa. The state of the collection is then always that which 
corresponds to the maximum number of complexions. 

Most of the calculations in statistical thermodynamics are based on 
Stirling’s approximation, which enables us to simplify the expression of the 
factorial logarithm n if the number n is large. It is written thus: 

In n ! = n In n - n = n In n [A.3 . 1 ] 


Modeling of Liquid Phases, First Edition. Michel Soustelle. 

© ISTE Ltd 2015. Published by ISTE Ltd and John Wiley & Sons, Inc. 
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A.3.1. The three statistical distributions 

Each element in the collection has an energy £}, and the number of 
elements which have that energy is The total number of elements will be 
N, such that: 

tf = 5>, [A.3.2] 


Hence, the total energy is: 

E = Y, n A [A.3.3] 

i 

The number of complexions, i.e. of configurations of the collection of 
elements, is represented by Q. 

The mean energy of an element is (e ) , and by applying relation [2.2], we 
find: 


(*) = f [A.3.4] 

Thus, for the number of objects in the state i, we find: 

n, = g, exp ( -a ) exp ( -fie . ) [A.3.5] 


where g i is the statistical weight or the coefficient of degenerescence or the 
multiplicity of the energy level £, . It is the number of different states with 
the same energy £ j . 


The coefficient j3 is a universal value which is: 


P = 


1 

kj 


[A.3.6] 


k B is Boltzman’s constant (the quotient of the joule constant, R, by 
Avogadro’s number N a ). 
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Depending on the nature of the molecule-object, three kinds of statistics 
have been developed. 


A.3.1.1. Maxwell-Boltzmann statistics 


Maxwell-Boltzmann statistics applies to objects for which there is no 
room to resort to quantum mechanics, i.e. fairly large and easily-discemible 
objects. This branch of statistics is also applicable to discernible or localized 
quantum objects, such as the “molecules” distributed at the nodes of a 
crystalline lattice. 

The coefficient crfrom relation [A.3.5] is given, in this case, by: 


exp (-o:) 


N 

2>/ ex pi-fa) 


[A.3.7] 


The distribution law becomes: 


Ngj exp (-/?£, ) 
exp (-/?£,) 

i 


[A.3.8] 


/? is defined by relation [A.3.6], 

A.3.1.2. Bose-Einstein quantum statistics 

Bose-Einstein quantum statistics applies to non-localized quantum 
objects, i.e. which are indiscernible and have integer spin (most of the 
molecules and the ions, the atoms). The distribution of the objects obeys the 
expression: 


g. exp (-or 
1-exp (-ar-/fe,) 


[A.3.9] 


The value of the coefficient cris difficult to determine. We shall come 
back to this point later on (relation [A.3. 1 1]). 
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A.3.1.3. Fermi-Dirac quantum statistics 

Fermi-Dirac quantum statistics applies for non-localized quantum 
objects, meaning indiscernible objects with fractional spin (some molecules 
and ions, and electrons). The distribution of the objects obeys the expression: 


g,.exp (-a- fie,) 
1 + exp (-a- /3e,) 


[A.3.10] 


The value of the coefficient a is as difficult to calculate as in the previous 
case. 

For the two branches of statistics for non-discemible objects, we content 
ourselves, for that coefficient a , with a limited expansion of the form: 


N a 

pO 

2 

,XY 

a 0 + a x 1 - a 2 


+ 


Z A 

v Z A j 


V Z A y 


+ ... 


[A.3.11] 


Laborious calculations show that the coefficients a; in this expansion are: 

[A.3.12] 


a o 0 , a x 1 , a 2 — 23/2 ’ fl 3 ^ ^ 3/2 ’ • • ' 


In the coefficient a 2 > the + sign applies for Fermi-Dirac statistics, and 
the - sign applies for Bose-Einstein statistics. 


A.3.1.4. Classic limit case 

The three quantical statistics (MB, BE and FD) come together to form 
one, called the classic limit case, if the following condition is fulfilled: 

exp (-or) «1 [A.3.13] 

Thus, or must be large. 

In these conditions, the three laws come together in the form: 


n,=g i exp(-a-/3e i ) 


[A.3.14] 



Appendix 3 211 


We can see that this limit case maintains the formula of the Maxwell- 
Boltzmann distribution. Thus, a is determined by the relation: 

exp(-a A ) = ^ [A.3.15] 

Z A 

We can show that in this case, we also have the condition: 

N«g t [A.3.16] 

i.e. if the number of particles is much less than the number of possible states. 


A.3.2. Partition functions of a molecule object 
A.3.2.1. Definition 

The partition function of a molecule object of a collection is the sum z 
defined by: 


f ^ 


z =Z^ ex P =Z^- ex p(“^) 


V k B r 1 i 


[A.3.17] 


The sum is extended to all the levels of energy that the object may attain. 


A.3.2.2. Independence of the energies 

The complete partition function for a system includes terms that refer to 
the various forms of energy: nuclear, electronic, vibration of the molecules, 
their rotation, their translation and the energy of interaction between the 
different molecules. 

To simplify, we accept that these different forms of energy, for a 
molecule, are independent. 

In these conditions, we can write the total energy of a molecule as 
the sum of the different contributions of the forms of energy - nuclear 
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£ n , electronic f , vibrational e v , rotational £, , translational £ t and 
interactional £, , so that: 


£-£„+ £ e + £ v + £,. + £,+£, 


[A.3.18] 


The molecule’s partition function becomes: 


z=Yj ex p 

i n 

Z ex p 


f £ : X 


V k B T ) i e 


f £ : ^ 


Z ex P X cx p 


V k B T JK 


\ k B r J 


( £, ^ 


V k B T j i, 


( £. ^ 


Z ex P -irV 2 >p 


v k B^y 


f £, x 


[A.3.19] 


v V; 


Thus, we reveal partial partition functions, which are afferent to the 
different forms of energy: 


=Z ex p 

i n 

z,=Z cx P 


f £, ^ 


v k s :r y 


f ^ 


v k B^y 


^,=Z cx p 

i e 

, z , =Z ex p 


f £, ^ 


V 


f £, ^ 


v k e^y 


>A = Z ex P 

i v 

> z i =Z ex p 


f £; ^ 


V k B T J 


f £: ^ 


v. k B^y 


[A.3.20] 


The overall partition function then assumes the form of a product of the 
partial partition functions: 

Z = Z n Z e Z V Z r Z t Z I [A.3.21] 

Sometimes, we use the hypemym internal contribution to speak of the 
product: 

Z mX ~ Z n Z e Z v Z r [A.3.22] 

It is the product of all the contributions other than those of translation and 
interaction. 

Thus, obviously, the overall molecular partition function becomes: 


[A.3.23] 
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A.3.2.3. Partial molecular partition functions, relating to the different 
motions 


By applying definition [A. 3. 20], we can calculate the contributions of 
each of the motions of the molecule to the molecular partition function. 

A.3.2.3. 1. Translation 

The molecule has three degrees of freedom for translation. We can show 
that if it is subject to no other constraint than having to remain in the 
volume V, the contribution of translation is: 


z, = r 


2nmkfT 


[A.3.24] 


For the perfect gas, with no interaction (z/ = 1 ) between the molecules, 
the molecular partition function can be written as follows, in light of 
relation [A.3.24]: 




2/nnk B T 


, 3/2 


[A.3.25] 


Thus, the translational partition function of the perfect gas is: 


Z Hpf)= V 


2nmkfT 


, 3/2 


[A.3.26] 


A.3.2.3. 2. Rotation with moment of inertia I 

A molecule may have two or three degrees of freedom for rotation. We 
distinguish three families of molecules. 

For homonuclear diatomic molecules (i.e. where the two atoms are 
identical), the partition function per degree of freedom is: 

4n 2 A B T 

— h 5 


z. 


[A.3.27] 
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For heteronuclear diatomic molecules, the partition function per degree of 
freedom is: 


z 


r 


4n 2 Ik B T 
— 1 ? 


[A.3.28] 


In the case of complex molecules, the partition function per degree of 
freedom is: 


8 n 2 (inkjf 2 VWs 
ah 2 


[A.3.29] 


where a is a coefficient of symmetry which depends on the complexity of 
the molecule, whose value is at most a few units. 


A.3.2.3.3. Vibration of frequency V 

We can show that the partition function of one degree of freedom of 
frequency vis given by: 


z 


V 


exp 


hv 

2k7 


1-exp 


hv 


[A.3.30] 


Note that if 


hv 

kj 


» 1 , we can content ourselves with a simpler form: 


z 


v 


Kf 

hv 


[A.3.31] 


A.3.3. Canonical partition function 

To use statistics for the characterization of the phases, it is interesting to 
construct the canonical ensemble. 
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A.3.3.1. Canonical ensemble 

A so-called canonical ensemble is an ensemble, comprising replicas of 
the system under examination. Each element is closed, so the number of 
molecules N is identical in every element of the ensemble. This number is 
constant because there is no exchange of matter between the elements, or 
between those elements and the outside world. The volume V is the same for 
all the elements. The elements are in thermal contact with one another, and 
can thus exchange energy. Their temperature T is identical. Each element has 
an energy Ej. The total energy of the canonical ensemble would be E c . This 
energy is constant because the ensemble is isolated from the outside. 

A.3.3.2. Canonical partition functions 

Similarly, as for molecules, we define the partition function for the 
canonical ensemble by the sum: 

Z c =£ exp (-/?£,) [A.3.32] 

j 

This sum is extended to all the elements in the ensemble. 


A.3.3.3. Canonical partition function and molecular partition functions 

We can link the canonical partition function, firstly to the molecular 
canonical functions, and secondly to the thermodynamic functions which 
define the phase on the macroscopic level ( U , F, G, S, etc.). These two types 
of relation use the canonical partition function to form the link between the 
microscopic definition of the phase and its macroscopic thermodynamic 
properties. 

In order to calculate the canonical partition function on the basis of the 
molecular functions, we distinguish two cases, depending on whether the 
molecules are discernible or indiscernible. 

A.3.3.3. 1. Case of collections of discernible molecules 

If the molecules are all identical and discernible, we show that we can use 
the expression: 

Z c =z N 


[A.3.33] 
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If the phase comprises several types of molecules - N A molecules of A, in 
N b molecules of B, etc., the canonical partition function becomes: 

Z c =z N f.z N f.... [A.3.34] 

A.3.3.3.2. Case of collections of indiscernible molecules 

If the molecules are all identical and indiscernible, we show that we can 
use the expression: 

N 

Z c =— [A.3.35] 

c N\ 

If the phase comprises several types of molecules - N\ molecules of A in 
N B molecules of B, etc. - the canonical partition function becomes: 

N A n b 

Z c = ^-.-5-... [A.3.36] 

L M I M I 
JV A • JV B- 

In the case that the medium under examination is a mixture of several 
phases a, fi y etc., the canonical partition function of the ensemble is the 
product of the canonical partition functions of the different phases, in 
accordance with: 

Z c =Z ( " ) .Z ( c p) .Z ( f ) .. [A.3.37] 

Each canonical partition function for each of the phases obeys one of the 
relations [A.3.33], [A.3.34], [A.3.35] or [A.3.36], 


A.3.4. Interactions between molecules 

Let us consider the mixture of N independent molecules, with no 
interaction between them. This collection is treated as an assembly of non- 
discemible quantum particles in the classic limit case of validity of relation 
[A.3.13], 

If we use the notation (Oj to denote the set of position coordinates (jti, x 2 , 
Xi)i of the particle i, between two particles i and j there is an energy of 
interaction £ t j, and we suppose that the interaction energy only contains 
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terms such as £ t j, although many molecules interact mutually. This 
hypothesis is known as the pairwise interaction model. The configuration 
integral thus takes the form: 

h = j rjn ex p-i% d < t A - 3 - 38 ] 

Nv f i<j Al' 

The integrals are extended to the available volume for the molecules, so 
Nvf. 


We suppose that there is a zone of interaction for the molecule, the sphere 
of influence, and that there is only one molecule j which is within the sphere 
of influence of a molecule i. Thus, by ignoring the volume of those 
molecules in relation to the total volume, we establish a term of interaction 
in the molecular partition function of the form: 


I, = exp- 


N 2 B aa (T) 

V 


[A.3.39] 


The term B 4A (T) being given by the expression, a function of the 
distance r between two molecules: 


B 


i (T) = - ] -\A7lr 


exp 


KT; 


-1 


dr 


[A.3.40] 


Usually, we content ourselves with a limited expansion of this 
contribution in the form: 


N 2 B aa (T) 

V 


[A.3.41] 


For a collection of A molecules of gas (indiscernible molecules), this 
gives us a canonical partition function which has the expression: 


ZlL 

A! 


1- 


N 2 B m {T) 

V 


f InmkJ ’ 

h 2 J 

3 / 2 " 

N 

N rrN 

A! 


1- 


N 2 BJT) 


[A.3.42] 


V 
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The second coefficient of the virial is therefore: 
B 2 {T) = N a B A4 (T) 


[A.3.43] 


The energy e AA of interaction between two molecules may take a variety 
of forms, particularly the Lennard-Jones form, which is due to the van der 
Waals forces: 





= — 2 


v d J 


+ 


( d ^ 
u 0 

d 


[A.3.44] 


where do is the distance between two molecules which corresponds to the 
minimum. 


Up until now, we have defined an energy of interaction between 
molecules £y. It is possible to work another way, considering the molecules 
without interaction and introducing a term representing the overall energy of 
interaction E h 


The configuration integral is written, at the level of the canonical 
partition function, by replacing the formula [A.3.38] by: 


4JJ 


exp 


k Jj 


d (o' 


[A.3.45] 


A.3.5. Canonical partition functions and thermodynamic functions 

We can demonstrate the following formulae to express the 
thermodynamic functions on the basis of the canonical partition function. 


For the internal energy: 


x — 1 9 In Z ct i ^ 9 In Z a a\ 

U(T) - U( 0) = -X .J fA) = k B T X ~ 

A Op A 


dT 


[A.3.46] 
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For the entropy: 


S = k B 


lnZ c 


1 3lnZ c 
T 3 ln/7 


[A.3.47] 


For the Helmholtz energy function: 

^lnZ C(A ) _ 

F(T) - F( 0) = = -k B ZX lnZ C(A) [A.3.48] 

P A 

NOTE (on the molar values).- In the expressions of the molecular and 
canonical partition functions, we worked for a certain number of molecules - 
a number of particles symbolized by N A for component A. 

To obtain the molar values of the thermodynamic functions, it is wise to 
choose Avogadro’s number (N a ) as the value of N A . 

To obtain the value of a function for an amount n A (in moles) of the 
component, then for N A we use the product n A N a . 


A.3.6. Equilibrium constants in the liquid phase and partition functions 


As the thermodynamic constants - particularly the Gibbs energies - can 
be expressed on the basis of the partition function, the same must be true of 
the equilibrium constants. 


The equilibrium constant of a reaction is expressed on the basis of the 
partition functions of the reagents and products, in the case of liquid perfect 
solutions of dilute ideal solutions, by the expression: 




A r h\T) 
R T 


[A.3.49] 


The term z° (i) is the molecular partition function of the component i, the 
partition function expressed in relation to a molecule, so: 


Z m(i) ~ 


N 


«N a 


[A.3.50] 
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The term A r h°(T) is the linear combination, weighted by the 
stoichiometric numbers, of the residual vibrational energies of each of the 
substances. 

Thus, if each component has kj vibrational degrees of freedom with the 
fundamental frequency v° k , we would have: 

A,. U °(0) = ^XkZ<) = M°(0) [A.3.51] 

i k , 

If the solution is not perfect, this term is joined by the energy of 
interaction of mixing at the temperature of 0 K. 
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